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Recently, a family of exact force-free electrodynamic (FFE) solutions was given by Brennan, Gralla 
and Jacobson, which generalizes earlier solutions by Michel, Menon and Dermer, and other authors. 

These solutions have been proposed as useful models for describing the outer magnetosphere of 
conducting stars. As with any exact analytical solution that aspires to describe actual physical 
systems, it is vitally important that the solution possess the necessary stability. In this paper, 
we show via fully nonlinear numerical simulations that the aforementioned FFE solutions, despite 
being highly special in their properties, are nonetheless stable under small perturbations. Through 
this study, we also introduce a three dimensional pseudospectral relativistic FFE code that achieves 
exponential convergence for smooth test cases, as well as two additional well posed FFE evolution 
systems in the appendix that have desirable mathematical properties. Furthermore, we provide an 
explicit analysis that demonstrates how propagation along degenerate principal null directions of 
the spacetime curvature tensor simplifies scattering, thereby providing an intuitive understanding of 
why these exact solutions are tractable; i.e. why they are not backscattered by spacetime curvature. 


I. INTRODUCTION 

Force-free electrodynamics [IHS] (FFE) is a simplifica¬ 
tion to the joint electromagnetic and plasma dynamics 
that is applicable in the limit of magnetic domination. 
Within FFE, the inertia of the plasma is neglected, so 
that the equation of motion for the plasma is not required 
to close the set of evolution equations. This leads to a sig¬ 
nificant reduction in computational complexity, with the 
presence of the plasma becoming a nonlinear modifica¬ 
tion to the vacuum Maxwell equations. Astrophysically, 
FFE is recognized as the appropriate limit for describing 
the magnetospheres of black holes [T] and particularly 
neutron stars [2], where one can find intense magnetic 
fields of 10® - 10^® gauss accompanied by charged parti¬ 
cles supplied by electron-positron pair-production mil. 
FFE is an integral part of most proposed mechanisms 
for extracting rotational energy from neutron stars [2] or 
black holes mm, and electromagnetic dominance (over 
gas dynamics) is argued to be valid in all ultra-relativistic 
outflows [6]. For instance, the jets in quasars and active 
galactic nuclei [7] or gamma-ray bursts [5] are generally 
simulated using the FFE approximation. To understand 
these astrophysical phenomena, it is therefore important 
to analytically (e.g. mEHn]) and numerically (e.g. [13- 
[75] ) study the solutions to the EEE equations. 

One step in this direction was achieved recently with 
the presentation of a family of analytical solutions in 
Kerr spacetime by Brennan, Gralla and Jacobson [IT] , 
which combines and generalizes some earlier solutions 
by Michel [101 177] . Menon and Dermer |7], and puts 
them in a language more accessible to relativists. It has 
been suggested that these solutions (containing nonlin¬ 
ear ingoing or outgoing waves) can describe astronomi¬ 
cal systems such as the outer magnetosphere of a pulsar 


[171 l30l ITT] , describing the mechanism for transporting 
energy extracted from the interior regions towards infin¬ 
ity. Because the EFE equations are highly nonlinear, an¬ 
alytical solutions are relatively rare (see [37U35] for some 
additional solutions in extremal Kerr spacetime), so it is 
worthwhile to examine these solutions in greater details, 
especially those aspects related to their applicability to 
real astronomical problems. For the benefit of finding 
further exact solutions, it is also interesting to study the 
properties that make these known solutions tractable. 

A remarkable feature of these solutions is that their 
wave contents are not backscattered^ by the spacetime 
curvature, a fact that significantly simplifies the analysis, 
and in no small part contributes to the possibility of ex¬ 
pressing these solutions in closed form. Such scattering- 
avoidance behaviour is not typical among waves trav¬ 
elling in a curved spacetime; generically the waves will 
scatter against the spacetime curvature and travel inside 
(as well as on) the null cones. In addition to not being 
scattered by spacetime curvature, these solutions are also 
not backscattered by nonlinear electromagnetic interac¬ 
tions, which makes them particularly efficient channels of 
energy transfer. A question then naturally arises: does 
the physical specialness and mathematical simplicity of 
these solutions equate to fragility? In other words, are 
these solutions stable under initial perturbations? The 
answer to this question is critical, as a negative answer 
would mean that these solutions do not describe realis- 


^ In this paper, we will refer to “backscattering” as the scatter¬ 
ing capable of altering the propagation direction of the waves, 
or introducing a Coulomb component where initially there was 


none. 
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tic astronomical systems that are always subject to per¬ 
turbing influences from their surrounding environments If 
these solutions represent “repellers” in the FFE solution 
space, then one would expect that a small perturbation 
in the initial data will quickly drive the physical system 
away from them and perhaps establish alternative, phys¬ 
ically less special and less efficient energy transportation 
channels. Even worse, if uncontrolled growth in the mag¬ 
nitude of the perturbations appears due to the presence 
of unstable modes, it would be a sign that the mathemat¬ 
ical model may be intrinsically inadequate for describing 
real physics. 

On the topic of stability of force-free and magneto- 
hydrodynamical conhgurations, there is a rich and im¬ 
pressive list of literature on the structure (see e.g. [36l - 
m]) and stability of jets. Researchers have examined in 
detail the various instability types that could alter the 
jet structure or even disrupt it. For example, Refs. gl- 
156] examined various types of current driven instabilities, 
providing important observations such as that instabili¬ 
ties can change the current density in the jets [sni[S2], 
forming structures that can heat and accelerate particles, 
and that increased magnetization tends to have a stabi¬ 
lizing effect m- Other instability types in jets such as 
Kelvin-Helmholtz and pressure driven modes have also 
been the subject of intense studies (see for example gH- 
|62]). Currently, a consensus on the reason behind the 
remarkable stability of observed and numerically simu¬ 
lated (see e.g. gnUBMO]) jets is still lacking in], and 
will continue to be a fascinating area of research. 

In this paper, we tackle a rather different problem. 
(1) We do not examine collimated jets, but rather more 
isotropic radiation which appears to contribute to a sig¬ 
nificant [26] or even dominant m portion of the energy 
budget in the radiation emitted by e.g. a binary black 
hole system inside a common magnetosphere. This type 
of radiation has received less attention previously, but the 
particular scatter-avoiding solutions we examine, if sta¬ 
ble (their stability has not been examined before), may 
prove to be a preferred (as it is the most efficient, with¬ 
out backscattering of energy flux) channel through which 
isotropic energy flux escapes through the magnetosphere, 
thus providing us with a perfect entry point for further 
research. Because a large pulse of isotropic radiation is 
emitted during the merger phase of a black hole binary 
[26] . which can potentially be picked up by observers on 
Earth, such research should have relevance for multimes¬ 
senger (gravitational and electromagnetic waves) astron¬ 
omy. (2) Because these solutions are envisaged to be the 
couriers that carry energy across magnetospheres of black 
holes or neutron stars, relativistic effects would become 
important, so our analysis will necessarily have to take 
spacetime curvature into account. In contrast, most of 
the previous studies on jet stability assume flat space- 
time, as jet disruption does not occur until very far from 
the central compact object. (3) We are examining global 
(while jet stability studies concentrate on the vicinity of 
the jets) solutions, and so there are subtle new instability 


types that may emerge. For example, there is no globally 
regular vacuum counterpart to the solutions we examine 
[HEa, so does that mean singularities similar to those 
seen in the vacuum case would generically develop even 
in the force-free case when we introduce perturbations? 
This subtle potential instability, whose underlying source 
is global in nature, would not have been included in the 
consideration of typical severe plasma instabilities. On 
the other hand, the solutions we examine do not have ex¬ 
treme features such as a jet boundary, so they are likely 
less prone to many instabilities that a jet would suffer. 
In fact, a cursory glance gave us no reason to strongly 
expect these typical plasma instabilities to severely im¬ 
pact these solutions (we caution however, due to different 
assumptions regarding the structure of the solution and 
the nature of the plasma, jet stability results may not be 
directly applicable to the present study). (4) Whereas 
previous studies mostly concentrate on the stability of 
a physical feature (i.e. jets), we examine the more re¬ 
strictive case of the stability of a particular family of 
exact analytical solutions. Therefore, even though phys¬ 
ical processes that slightly alter the exact structure of 
the jets - but do not disrupt it - may not be regarded 
as serious instabilities, they would in our case be neces¬ 
sarily recognized as a problem, as they would nudge the 
actual physical configurations away from being precisely 
the same as the exact solutions. 

A consideration of all possible instabilities to an ana¬ 
lytical solution is exemplified by the study of Kerr metric 
stability (see e.g. Ref. [73| for a summary). The prob¬ 
lem can be attacked by first studying the mode stability 
(see e.g. [SD] for an example in the context of jet stabil¬ 
ity studies) through solving linearized perturbation equa¬ 
tions assuming separable solutions. For example, a recent 
paper m studying perturbations of magnetic monopole 
and Blandford-Znajek solutions showed no unstable indi¬ 
vidual FFE modes^. It is then an arduous task to further 
prove linear stability, as it is not guaranteed that all lin¬ 
ear perturbations can be decomposed into such modes, or 
that the sum of inhnitely many stable individual modes 
remains stable m- A proof for full nonlinear stability 
is more difficult still. On the other hand, while a rigor¬ 
ous proof is currently out of reach, important evidence of 
nonlinear stability can often be found using fully nonlin¬ 
ear numerical simulations. For example, several studies 
of the numerical robustness of the Blandford-Znajek pro¬ 
cess can be found in Refs. HSHSMOj. Indeed, the ability 
to examine stability is cited as one of the principal mo¬ 
tivations for developing fully time-dependent numerical 


^ A similar analysis for perturbations over the FFE solutions we 
are interested in would be much more difficult. These solutions 
have a null current that couples to the perturbing fields just like 
the Blandford-Znajek case. However, unlike the current for the 
Blandford-Znajek solution, this null current is not proportional 
to a small parameter like the black hole spin, and cannot be 
treated using the perturbative techniques of Ref. ED. 
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FFE codes [T71 [7S] . In this paper, we adopt this more 
accessible numerical approach to studying nonlinear sta¬ 
bility. 

This paper is organized as follows: We begin by intro¬ 
ducing the force-free equations in curved spacetimes, as 
well as the details of the exact analytical solutions that 
we examine in Sec. m We then introduce in Sec. Ella 
new pseudospectral numerical code, used for our stabil¬ 
ity study. We also present several nontrivial tests demon¬ 
strating that the code achieves exponential convergence. 
Then in Sec. |IV[ we evolve a constraint-free perturbed 
initial data set, and show that the exact analytical solu¬ 
tions are in fact stable, despite their physical specialness. 
In Sec. [Vj we provide some derivations and arguments 
that provide an intuitive explanation as to why the ana¬ 
lytical solutions are not backscattered by spacetime cur¬ 
vature (this section is not necessary for understanding 
the numerics described in the earlier sections, and uti¬ 
lizes spinors extensively; therefore, readers who are only 
interested in the numerical aspects of this paper do not 
need to review Sec.|^. Finally, although not used in the 
numerical studies in this work, we also provide in Ap¬ 
pendix 1^ a couple of force-free evolution systems that 
have improved well posedness properties. 

In this paper, we adopt geometrized units with G = 
c = I, and use (—h -I—1-) for the metric signature. The 
beginning of the lowercase latin alphabet will be used 
to denote spacetime indices, and the middle of the al¬ 
phabet denotes spatial indices. Capital latin letters will 
denote spinor indices, while greek letters will index dif¬ 
ferent quantities in different sections, whose meaning will 
be clear from their context. Bold-faced letters will denote 
vectors and tensors. The numerical work in this paper is 
carried out within the pseudospectral code infrastructure 
of the Spectral Einstein Code (SpEC) [75] . 


II. THE FORCE-FREE EQUATIONS AND 
THEIR EXACT SOLUTIONS 

A. Some useful definitions 

In this paper, we will use the 3-1-1 form of the metric 

= —N'^dt^ + gij{dx^ + fi''dt){dx^ + fi^dt ), (1) 

where N is the lapse, (3 is the shift, and g is the (spatial) 
metric for the spatial hypersurfaces of constant t. The 
extrinsic curvature K of these spatial hypersurfaces is 
given by 


{dt-Cf,)g = -2NK, (2) 

which is a spatial tensor depending on both the geom¬ 
etry (metric) of the overall four-dimensional spacetime 
and the way we slice it. For example, the extrinsic cur¬ 
vature of a Minkowski spacetime can be nonvanishing if 
one picks unusual slicings. The operator on the left-hand 


side of Eq. ([^ is the derivative along the normal vector 
lt“ to the spatial hypersurfaces. 

When there is an electromagnetic field represented by 
the Faraday tensor F, we can also break it down into 
a 3-1-1 form, which will then allow us to write the force- 
free evolution and constraint equations in terms of spatial 
tensors in the next section. These equations will then re¬ 
semble those used in the numerical study of the Einstein 
equation, and can be handled with the same code infras¬ 
tructure. We define the electric and magnetic vectors 
as 

Fable ■ (3) 

Note that although we have written them as four-vectors 
in the definitions, they are really only spatial vectors with 
and vanishing. We will denote them as 3-vectors 
E and B below, with the understanding that a projec¬ 
tion into the spatial slices has been taken. Within the 
spatial slices, we will also use traditional vector calculus 
notations to simplify expressions, with, for example 

E-B = E^B^g,„ (ExB)* = /E^B'=e,,fc. (4) 


B. The evolution equations 


In this section, we write down a set of FFE equations 
with constraint damping capabilities that are numerically 
robust, although they possess some mathematically un¬ 
desirable properties that do not appear to hinder their 
performance in practice. The numerical studies carried 
out in the main body of this paper use this evolution sys¬ 
tem. Aside from this system of equations, we also provide 
in Appendix [^ two additional sets of equations with de¬ 
sirable mathematical properties, but which are numeri¬ 
cally less forgiving. We go through the derivation of these 
equations in some details for pedagogical reasons, as ex¬ 
isting literature tends to be brief and sometimes leaves 
out terms that should be included for curved spacetimes. 

We begin by writing down the Maxwell equations in 
curved spacetime, which are 


{dt - £p)E 
[dt - Cp)B 

V • E 

V B 


TV ATE -I- V X (iVB) - 47rAJ , (5) 

NKB - V X (NB ), (6) 

47rp, (7) 

0 . ( 8 ) 


To derive the current J, which is the spatial part of a 
four-current we begin with the force-free condition 
^o& j^4) _ states that the four-force density de¬ 

scribing the transfer of energy and momentum between 
the electromagnetic fields and the charged plasma parti¬ 
cles vanishes. This ensures that the stress-energy of 
the electromagnetic field remains dominant over that of 
the plasma. Indeed, we can derive the force-free condi¬ 
tion starting from the differential conservation of energy 
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and momentum VaT°‘^ = 0. Then when is the dom¬ 
inant contribution to we have mm 

- ^aT^M = = 0 . (9) 

In a 3 -b 1 decomposition, this translates into 

E-J = 0, pE-bJxB = 0, (10) 


where the second equation is the vanishing of the Lorentz 
force. To derive the force-free current J, we take the cross 


product between Eq. (10) and B which gives us 


B 

dTrA'^J = 47r7V(B • J) —- + AttN p 


E X B 

^2 


( 11 ) 


where the second term above can be seen as the charge 
density moving at the plasma drift velocity, and we can 
replace p with V • E/47r using one of the Maxwell con¬ 
straint equations. To further work out the current B • J 
along the B field, we note that 


= 0 ^ E-B = 0 or p = Q ,(12) 


and for nonvacuum solutions (vacuum here refers to 
= 0, with these solutions satisfying the force-free 
condition trivially) we would like to enforce the E • B = 0 
condition, which should be preserved along the timelike 
normal to the spatial hypersurfaces, and so 


{dt-Cp)^-^ = Q. (13) 


Using the definition of the extrinsic curvature tensor 
Eq. (|^, and substituting in the Maxwell equations, we 
obtain an equation for B • J that reads 


47r7VB • J = -E • V X (iVE) -b B • V x (fVB) 

-2NKijE^B^ + 2NKE • B , (14) 

(note that the extrinsic curvature terms on the second 
line appear to be missing in some of the existing liter¬ 
ature). Substituting Eq. (14) into Eq. 0 yields J in 
terms of E and B. Substituting this expression back into 
the Maxwell-equations yields the desired minimalist FEE 
evolution system. 

There is also a set of constraints that needs to be satis¬ 
fied, which comes from both the Maxwell equations and 
the force-free condition. When deriving the current J, 
we have explicitly used q — V • E/47r as the definition of 
charge density, so there is no need to enforce this con¬ 
straint. The nontrivial constraints are V • B = 0 and the 
force-free constraint E • B = 0. These two constraints 
are preserved automatically by the evolution equations 
m- Specifically, the V • B = 0 constraint is preserved 
by the original Maxwell equations and inherited by the 
force-free specialization. The E • B = 0 condition is 
also preserved, as we have explicitly used the condition 
(9t-/:^)E.B = 0 to derive the current. Physically, this 
condition fixes the magnitude of the conduction current 
along the B direction |T7], which would short out the E 


field along B by redistributing charge to eliminate the 
potential difference associated with that E component, 
thus enforcing E • B = 0. 

Although the V • B = 0 and E • B = 0 constraints 
are preserved by the evolution equations when they were 
satisfied initially, numerical noise inevitably creates some 
seed constraint violation that may grow further under 
the minimal evolution system. It is therefore benehcial 
to modify the evolution equations so as to be able to 
clean up the constraint violations as they emerge. For 
the E • B = 0 constraint, we adopt a strategy similar in 
form to Ref. m- Specifically, we add a damping term 
—A^i5(E • B)B/i3^ to dtE, so that the full set of evolution 
equations becomes 

|~^ 

{dt - Cp)E = NKE -b V X (A^B)- ■ E 

JD^ 

B 

-N— (B-VxB-E-VxE 

B2 ^ 

-2K,jE^B^ + 2KE • B -b (5E • B) (15) 
[dt - Cf3)B = NKB - V X (A^E). (16) 

The damping term is proportional to the constraint, and 
will not affect the physical constraint-satisfying solutions. 
However, it modifies the constraint evolution equation to 
a damped form 


(dt - Ci3)B ■ B = -SNE ■ B . 


(17) 


We note that our damping strategy differs from that of 
Ref. [27] , in that we have kept the original current terms 
from Eq. (14) and treated the new damping term as an 


addition instead of a replacement. In contrast. Ref. m 
removed all of the original current terms from the evo¬ 
lution equations, replacing them with only the damping 
term. With their strategy, the evolution equations be¬ 
come simpler, but those current terms forcibly removed 
from the evolution equations will resurface in the con¬ 
straint evolution equation, specihcally on the right-hand 
side of Eq. ( [l7| ). Therefore Eq. ( |T7| won’t reduce to 
{dt — C/ 3 )E • B = 0 when E • B = 0, so that constraint- 
satisfying FFE solutions at some instance of time can not 
stay constraint-satisfying as the simulation progresses. 
However, this does not invalidate the results of Ref. 
as the size of the constraint violation would be negatively 
correlated with S (and positively correlated with the mag¬ 
nitudes of the derivatives of B and E), as the damping 
will be activated when E • B grows too large. Indeed, 
Ref. |57| adopted a S greater than the inverse of the time 
step size and observed a well controlled E • B. Such a 
strategy however leads to a stiff term in the evolution 
equations that has to be treated with an implicit-explicit 
(IMEX) evolution scheme (37). As the analytical solu¬ 
tions we examine are exactly constraint-satisfying at all 
times, and since we do not use an IMEX scheme, we will 
use the damping term as an addition with a moderate 
coefficient S = 100, allowing us to enjoy its constraint 
cleaning benehts but avoid the aforementioned compli¬ 
cations. Finally, note that this damping term does not 
introduce magnetic monopoles, see Eq. (B51). 
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We note that the introduction of the additional damp¬ 
ing term replaces another alternative constraint cleaning 
strategy of removing the component of E along B after 
taking each time step, which has been widely utilized 
(e.g. Refs. [13 [m [T71 [77]). Such an alteration of 
the evolving fields at a discrete set of times (which de¬ 
pend on resolution) will result in a failure of the system 
to achieve the usual convergence behaviour that would 
be expected if the scheme were just applied to a set of 
differential equations without this alteration m- In con¬ 
trast, a damping term is less intrusive and its properties 
are more easily understood. We will show that our set 
of evolution equations displays the expected convergence 
behaviour in Sec. C 2[ In addition, because this damp¬ 
ing term does not contain derivatives, it will not affect 
the characteristic structure of the evolution equations, 
which is particularly helpful for our pseudospectral im¬ 
plementation. 

Lastly, we note that aside from the aforementioned 
E-B = 0 and V-B = 0 constraints, we have an additional 
constraint of < B^. When a region devel¬ 

ops, the plasma particles have to move faster than the 
speed of light in order to experience a vanishing Lorentz 
force. One can see this from the second term in Eq. 0 , 
which can be written as ^Vd, where Vd is the drift veloc¬ 
ity for the advection of the charge density [HIITj. The 
inequality > B^ then implies a super luminal Vd. For 
an even simpler demonstration, if we consider the special- 
relativistic point particle case, then the requirement for 
the vanishing of the Lorentz force ( 7 (E-fv x B) would im¬ 
ply |v| > 1 when > B^. Such super luminal motion is 
unrealistic, but FEE evolution equations cannot prevent 
it because the current (Eq. |11[ ) is derived without invok¬ 
ing plasma physics. Consequently, the FEE equations 
have no internal checks that enforce the B^ — E^ >0 
constraint. In other words, the B^ — E^ > 0 constraint 
is not strictly a constraint in the mathematical sense like 
V • B = 0 and E-B = 0, and solutions satisfying the 
force-free condition © are not automatically magneti¬ 
cally dominated. One symptom of the breakdown of this 
physical constraint is that, when strong waves interact 
m, the Alfven mode characteristic speeds become com¬ 
plex, breaking the hyperbolicity of the FEE evolution 
equations (because these are physical modes, augmenta¬ 
tions to the evolution equations like those in Appendix [B] 
will not be able to change their characteristic speeds and 
will thus not cure this hyperbolicity breakdown). An¬ 
other example is the formation of current sheets, in which 
the magnetic fields can vanish and then reverse direction 
m- Physically, the force-free assumption is invalid when 
B^ — E^ < 0. The actual plasma particles, which have 
inertia, will experience a nonvanishing Lorentz force and 
be accelerated. In other words, the system becomes dis¬ 
sipative averting divergences and possibly restoring 
magnetic dominance [13. A proper treatment of such 
regions would require special codes for the plasma [13 
that do not assume the force-free condition being met. 
For the numerical studies carried out in the later parts 


of this paper, we do not need such a sophisticated treat¬ 
ment as we are guaranteed magnetic dominance by the 
presence of a magnetic monopole in the solutions we ex¬ 
amine. We will discuss this in more details in the next 
section. In particular, we do not need to adopt the com¬ 
mon procedure of scaling down the E field after each time 
step to avoid electric dominance [I3II3II71|77|. 

We have now a basic set of evolution equations and 
the associated constraints. However, this does not auto¬ 
matically mean that we can simply plug them into the 
computer and run simulations, because we need to en¬ 
sure that the evolution equations are well posed. This 
is a rather technical discussion, which we have relegated 
to the appendixes. It turns out that the simple evo¬ 
lution systems given by Eqs. (15) and (16) (which we 
adopt for our numerical code) is not strictly strongly hy¬ 
perbolic (see Appendix but the violation is insignifi¬ 
cant enough that it does not create problems in practice 
for the simulations carried out in this work. Neverthe¬ 
less, we have formulated two additional evolution systems 
that are not only strongly hyperbolic, but also have ad¬ 
ditional desirable properties in terms of well-posedness. 
One of them is strongly hyperbolic even when E • B 7 ^ 0, 
while the other has a particularly simple constraint evo¬ 
lution behaviour. These systems are given in Appendix 
[b| Unfortunately, they (also including a similar evolution 
system from Ref. [Z3) introduce a term containing the 
second derivative of B (see Eq. B51) into the evolution 
equation for V-B, making it sensitive to high-frequency 
noises (see Sec. B4|. In other words, these mathemati¬ 
cally more satisfying systems are numerically less forgiv¬ 
ing. Therefore, we take a pragmatic approach, and use 
the simple system as given by Eqs. (15)-(16) for the nu¬ 
merical studies presented in the main text. However, the 
mathematically improved systems may yet prove useful 
for application in finite-difference and finite-volume FEE 
codes, which tend to be more forgiving in the presence of 
high-frequency noise than pseudospectral schemes. 


C. The exact analytical solutions 

In this section, we introduce the analytical solutions 
whose stability properties we seek to examine. The solu¬ 
tions are introduced in Ref. [13, and readers interested 
in their derivation should consult that reference. 

We begin by emphasizing that the spacetime back¬ 
ground used in this paper is the Schwarzschild spacetime, 
and the situation in a spinning black hole case is beyond 
the scope of our current study. This restriction stems 
from the need to enforce the magnetic dominance condi¬ 
tion that we discussed in Sec. El The scatter-avoiding 
wave-only analytical solutions given in Ref. m are the 
so-called null solutions, which can be seen as generaliza¬ 
tions to plane waves (see later in the section), and as 
such have E^ = B^ just like the plane waves. When we 
add a small perturbation to it, we can easily end up with 
E^ ^ ^2 j-ggjQjjs -which can not be properly accounted 
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I B_null I / I B_mono I 



FIG. 1: The ratio |Bnuii|/|I3mono| for the null"*" solution stud¬ 
ied, where Bnuii and Bmono are the null and monopole con¬ 
tributions in Eq. ( |33[ ), while the norm is defined as |Bnuii| = 
VBnuii • Bnuii. Shown is a cross-section in the a;z-plane, which 
extends radially from i?_ = 1.95M to = 195M. 


for under the force-free approximation. However, this is 
not a problem because instead of these null solutions, we 
in fact study the stability of what we call null'*" solutions, 
where a magnetic monopole is superimposed onto the null 
solutions (this prescription is also provided by Ref. [12], 
and |3I|), so all the scatter-avoiding properties of the null 
solutions are simply inherited by the wave components of 
the null’'" solution, yet magnetic dominance is maintained 
even with the presence of perturbations (Figs.[^an d[|de- 
pict the basic structure of field quantities for the null"'" 
solutions, and we will discuss them in more details later). 
Subsequently, we restrict ourselves to the Schwarzschild 
spacetimes because only when the black hole is not spin¬ 
ning would such a monopole addition not interfere with 
the null wave component, and the total solution (as a 
simple superposition of the wave and monopole compo¬ 
nents) still satisfy the force-free conditions. On the other 
hand, these two components would couple in a nontrivial 
way in a spinning black hole spacetime, and we do not 
have a simple null"'’ solution in that situation to study 
(the null"^ solution in that case has not been worked out 
yet). 

We note that because in astrophysical applications, we 
would have a split monopole background present in addi¬ 
tion to the energy-carrying waves, we are more interested 
in the stability of null"'" solutions instead of the null ones 
on physical grounds. So it is the stability of the null"'" 
solutions and not that of the null solutions that is the 
intended target of study in this paper. The study of the 
stability of null"'" solutions is not a surrogate for that of 
the null solutions. In fact, our stability results would 
unlikely translate from the null"'" to the null solutions, 
as we will discuss later in the Conclusion section. Fur¬ 
thermore, for our numerical study, we will pick the null"'" 



FIG. 2: Top left figure shows the electric field lines, and the 
top right figure shows the magnetic field lines for the same 
nulH solution as shown in Fig. the lines are coloured by 
and respectively. The bottom figures plot the time 
variation of the component, showing the wave propagating 
inwards. 


solutions with an ingoing (instead of outgoing ^) wave 
component, as these are regular on the future horizon 
of the Schwarzschild black hole. The outgoing solutions 
face the same possible “fragility due to specialness” issue, 
but need to be glued onto an astronomically realistic in¬ 
terior solution (not currently available); otherwise, they 
will represent energy flux emerging from the past horizon 

El- 

Before we dive into the details of the null+ solutions, 
we first introduce some formalism that will be needed. 
The structure of the solutions is most explicit in the 
Newman-Penrose (NP) formalism, wherein we express 
tensors under the NP null tetrad {l,n, m, m} that con¬ 
sists of two real null vectors 1 and n usually chosen to be 
in the outgoing and incoming directions, and two com¬ 
plex null vectors m and ni. This tetrad can be seen as 
the null version of an orthonormal tetrad m, relating to 
it via a rigid transformation. Under the null tetrad, the 
metric is a constant matrix 


/ 0 -1 0 0 \ 

I -1 0 0 0 I 
0 0 0 1 ’ 
Vo 0 10/ 


(18) 


^ The outgoing waves are obtainable through the transformations 
in Sec. 5 of Ref. |12| . 
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just as the metric is a constant matrix diag{ —1,1,1,1} 
under an orthonormal tetrad. The freedom for choosing 
the null tetrad is then given by the transformations that 
preserve the metric, namely the Lorentz transformations. 
Such a tetrad is just a basis to express components of a 
tensor in. In particular, the components of the Faraday 
tensor F under the Newman-Penrose null tetrad are the 
Newman-Penrose scalars /pQ, pi and p 2 , defined by 


00 

= FabFm\ 

(19) 

Pi 

= ^Fab{Bn^ + rh°^m^) , 

(20) 

p2 

= Fab'ffFrP'. 

(21) 


the Faraday tensor as 

Fab = 43 ? {PoFi[anb]) . ( 22 ) 

which can be shown to satisfy 

^F^^Fab = -^*F^^*Fab = B^-E^ = 0, (23) 

^Fab*F^^ = E-B ==0. (24) 

These conditions form the definition for a null wave. We 
have the following observations for our ingoing null wave 
solutions (which are also inherited by the wave compo¬ 
nent of the null’*' solutions) 


These three complex numbers are simply the 6 real com¬ 
ponents of Fab (three from E and three from B, and these 
vectors are of course just F’s components written in the 
coordinate basis) recombined. 

One reason why the NP tetrad is nicer than a coor¬ 
dinate basis is that we can pick the 1 and n in it to be 
pointing along special directions such as the outgoing and 
ingoing null directions of a Schwarzschild spacetime (we 
will be using the so-called Kinnersley tetrad [HI], which 
has this property. See Appendix for their detailed 
expressions). This specialness in the basis orientation 
translates into special meanings for the components of 
Fab associated with these bases, and subsequently po, Pi 
and p 2 can be interpreted as the incoming wave. Coulomb 
background and outgoing wave pieces of F, respectively. 
In particular, this identihcation is physical and gauge (co¬ 
ordinate choice) invariant, as the outgoing and ingoing 
directions as identihed by 1 and n in a Kinnersley tetrad 
are those determined by the spacetime geometry and are 
unambiguous [H0|j and not based on the radial direction 
of arbitrary coordinate systems, which can be subject to 
random coordinate transformations that change the ra¬ 
dial direction. 


As E and B and the three NP scalars are simply F com¬ 
ponents written in different basis, it is not surprising that 
we can translate the force-free equations (15) and (16) 
into equations for p^^ pi and p 2 , and that the equations’ 
structure clarifies tremendously because all the quanti¬ 
ties in it now have clear physical meaning. This is the 
approach taken in |12j . 

For the ingoing null+ solutions that we are interested 
in for this work, we need outgoing p 2 = 0, while the 
ingoing wave pQ and background pi are nonvanishing. It 
turns out however that in a Schwarzschild background, 
the equations for po and pi decouple, so we can solve for 
po first, assuming pi = 0, and then add a monopolar pi 
back later. The ansatz Ref. m uses when solving for pQ 
is then pi = 0 = p 2 , s^s well as an additional condition 
(guess) that the current flows in the ingoing null direction 
n. Under these assumptions, we get an equation for po 
we can actually solve. 

After obtaining an expression for po (like the one we 
will write down later in the section), we can reconstruct 


1. Locally, the two null wave conditions imply that the 
E and B fields are orthogonal to each other, and 
share the same amplitude just like a plane wave. 
These solutions are of course not really plane waves, 
see the E field lines in Fig. (also note the B field 
lines in that figure are for null+ and not null solu¬ 
tions, and are not applicable for our present discus¬ 
sion). 

2. The amplitude of E and B, and thus the wave, is 
given by po- 

3. This “local plane wave” travels along the ingoing 
n direction (radially inwards), which is also the di¬ 
rection of the current. 

4. The presence of the current is an extra freedom not 
available to vacuum Maxwell equations, so while 
FFE null solutions can be globally regular, singu¬ 
larities would develop in a vacuum null solution 

m- 

5. pi and p 2 vanish everywhere in this solution, mean¬ 
ing that the null wave is not backscattered (in 
which case an outgoing component would be gen¬ 
erated) by either the current or the spacetime cur¬ 
vature. And this property is shared by outgoing 
solutions, which makes them highly efficient chan¬ 
nels for transferring energy. 

6. The magnitude of po (and thus the amplitudes of E 
and B) is not limited. The solution represents fully 
nonlinear waves satisfying the nonlinear FFE equa¬ 
tions. This allows the outgoing variant of these so¬ 
lutions to carry an arbitrarily large energy flux from 
the interior region through the magnetosphere. If 
energy flux travels via scattered solution on the 
other hand, we would likely see much of the en¬ 
ergy turn back and be swallowed by the black hole. 
This is not what simulation of binary merger in 
magnetospheres appears to show. 

7. When we add a monopole background, and take 
the small amplitude/linearized limit, these waves 
reduce to the travelling waves discussed in Ref. [H2 . 
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And if we further take the eikonal limit such that 
the wavelength is much smaller than the radius of 
curvature, these waves become Alfven waves trav¬ 
elling along the background magnetic field [55] , 

Now that we have a broad picture of what these so¬ 
lutions look like, we turn to the details of a particular 
representative example that we will study in this work. 
There are infinitely many solutions to the equation, 
and we cannot numerically examine the stability of each 
and every one of them. So instead, we pick an arbi¬ 
trary representative solution which does not have any 
special properties (it shares the same basic current/field 
structure and other physical properties listed above with 
all of its siblings, and possesses no special symmetries) 
that would make it more stable than any other solution. 
Therefore its stability should be seen as strong evidence 
that most, if not all, of the other solutions obtained in 
the same manner should also be stable. 

One such representative solution is give in Sec 4.2.4 of 
Ref. [H], whose is given by 


= 


Ap 


(25) 


where 


A = — 2Mr and p = — (26) 

r 

in the ingoing Kerr (Eddington-Finklestein) coordinates 
with M being the mass of the background 
Schwarzschild black hole. The quantity is a real func¬ 
tion of the form 


/^ = 15F(z/) sin^ 0cosdCOS')/, (27) 

where F(v) is an arbitrary function specifying essentially 
the time dependence of the solution (y is the null coordi¬ 
nate in ingoing Kerr). The quantity also determines 
a companion function 

f = —J (d^ndo. 

smO J 

= —5F{i/) sin ij} sin^ 6 . (28) 


Given these expression we have now 


sin^ 6 


4>o = — [5F{i/) {3 cos 9 cos tjj — i sin'll))] , (29) 


while the current is 
1 


J = 


27ry2A 


\2QF{v) COS'!/sin0(2cos^ 9 — sin^ 6*)](>)0) 


flowing along the ingoing congruence tangential to the 
ingoing null base n. 

When we numerically implement this solution, we will 
use the Kerr-Schild coordinate system {t,x,y,z), in¬ 
stead of the ingoing-Kerr coordinates originally utilized 


in Ref. [T2]. In Kerr-Schild, we have 
3(a;^ -I- y'^)F{r 1) 


4>o — 


Szsin^tan ^ 
1 


- -7rsgn(a;) 


+ir cos ftan ^ ^ — -7rsgn(a:) 


(31) 


where r = +y^ + We can also now add a mag¬ 

netic monopole piece 

= (32) 

2 

(it has the same expression in ingoing Kerr and Kerr- 
Schild coordinates) to obtain our final magnetically dom¬ 
inated null+ solution. We can do this (despite the equa¬ 
tions being nonlinear) because the null and monopole 
solutions decouple in a Schwarzschild spacetime in the 
sense that the monopole has no currents, and its field 
tensor F does not exert a force on the current of the null 
solution [3T] . Therefore, the superposition of the null and 
monopole solutions still satisfy the force-free condition, 
and yield another FFE solution. (We note that although 
a magnetic monopole is not physically realistic, a more 
realistic split-monopole solution can be constructed by 
gluing two copies of the nulR solutions together, with a 
current sheet on the interface [301 El].) The monopole 
only contributes to the Coulomb part, and the field ten¬ 
sor is now given by 

Fab = 43? {(l)ofh[anb] + (/irufaTObj) . (33) 

When combined with the expressions for the tetrad basis 
summarized in Appendix [^ one can generate the Fara¬ 
day tensor and subsequently the E and B fields via 

Ea = F^^Tb , = (l/2)e“'’‘=‘'FcdTd (34) 

where T is the timelike normal to the Kerr-Schild spatial 
slices. We do not reproduce the full expressions for E 
and B in Kerr-Schild coordinates, as they are long and 
tedious. We can however, plot figures that further illus¬ 
trate the properties of the null"^ solutions. 

For concreteness, we first need to pick the parameters 
in our (/o and (j)\. Because (j)i drops off as ^ 1/r-^ (see 
Eq. [Mt while (/o drops at a slower rate of ~ 1/r (see 
Eq. |29p , we pick q = 1000 and some arbitrary values (so 
the solution has no special stability properties as com¬ 
pared to the rest in the family) F{v) = Acos(n^) with 
A = 1 and = 0.1 (i.e. the solution is time-dependent), 
so that the monopole is large enough to ensure mag¬ 
netic dominance at the outer edge of our computational 
domain (see the next section). The ratio between the 
null and monopole contributions to B is shown in Fig. [l] 
where we see that the two contributions are comparable 
in magnitude in the outer regions of the computational 
domain. 

In Fig. we plot the field lines for E and B for the 
nulR solution with the parameters set in the last para¬ 
graph. The top left panel of the figure shows the E field 
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lines, which clearly demonstrate that the solution is fully 
three dimensional, without any axisymmetry. The top 
right figure shows the magnetic field lines, which is clearly 
affected by the background monopolar contribution. The 
bottom panels show the time variation of the com¬ 
ponent of the electric field, demonstrating the time de¬ 
pendence of our solution and the fact that the wave is 
propagating inwards. Looking back at Fig. one also 
sees that the wave component of this particular null'*" so¬ 
lution is rather isotropic and not concentrated near the 
poles so there are no jet-like features. 


III. A PSEUDOSPECTRAL NUMERICAL CODE 


In this section, we briefly introduce the pseudospec- 
tral code used for evolving the force-free equations. As 
is noted in Refs. HZ], a pseudospectral code is espe¬ 
cially suited to the task of examining stabilities, as it 
avoids erroneous instabilities that may be triggered by 
the larger numerical noise present in finite-difference or 
finite-volume schemes. An example is pointed out by 
Ref. |T7j , that during a study of the Sweet-Parker recon¬ 
nection in Ref. [13], a spectral code was observed to not 
suffer from secondary island formations resulting from a 
tearing-mode instability, in contrast to results obtained 
using finite-difference and finite-volume schemes. 

In addition to our main evolution code, we also imple¬ 
ment an initial data solver to ensure that the constraints 
are properly satisfied. We note that this is an additional 
improvement relative to existing codes, as many previous 
studies using finite-difference or finite-volume codes have 
evolved initial data that violated the constraints (they 
however employ constraint cleaning schemes to remove 
the violation later during evolution). 


A. The code infrastructure 

In the interest of completeness, we briefly introduce the 
basic infrastructure employed in this work. Our force-free 
code is a module of the SpEC code, which was developed 
primarily to study binary black holes in full general rel¬ 
ativity. Readers interested in the details are encouraged 
to consult m and the list of articles shown there. 

The basic setup of the computational domain used in 
this work is a spherical shell (see Fig. [^, whose inner 
edge is at i?_ = 1.95M, and whose outer edge extends to 
i?_i_ = 195M, where M is the mass of the Schwarzschild 
black hole, which sets the length and time scales for the 
simulations, and will be used as a unit in the plots. The 
inner edge terminates just inside the event horizon at 
2M (the semitransparent surface in Fig. [^, effectively 
excising the singularity from the computational domain, 
so we won’t run into related numerical issues. The inner 
edge being inside the event horizon means there will be 
no information coming through that boundary into the 



FIG. 3: Half of the computational domain and the spectral 
grid (spectral collocation points are at the intersection of the 
black lines), the semitransparent sphere represents the event 
horizon of the Schwarzschild black hole. 


computational domain, and so we do not need to impose 
boundary conditions there. 

For the outer boundary, the evolution variables E and 
B are first translated into the characteristic modes (see 
Appendix for their detailed expressions), which can 
be seen as waves propagating normal to the boundary, 
carrying information into and out of the computational 
domain. To ensure there are no additional incoming per¬ 
turbations during the simulation, which would create the 
false impression of instabilities, we use the analytical ex¬ 
pressions for E and B from the exact null'*' solutions 
to compute the clean incoming characteristic modes and 
impose them as boundary conditions (so there is no ad¬ 
ditional perturbation being carried by these waves into 
the computational domain). On the other hand, there 
are no conditions imposed on the outgoing characteristic 
modes. This way, there is no new perturbations coming 
into the computational domain, but the perturbations 
already inside are allowed to exit. 

In order to carry out parallel computation, we break 
the entire computational domain into eight concentric 
spheres, whose boundaries are seen as dense concentra¬ 
tions of black circles in Fig.[^ Each subdomain is handed 
to a separate processor for computations. The commu¬ 
nication between subdomains is also done via character¬ 
istic modes, where the outgoing characteristic modes of 
one subdomain are matched onto the ingoing character¬ 
istic modes of its neighbouring subdomains via a penalty 
method [55H55] so that any discontinuity across the sub- 
domain boundaries is forced to vanish over time. Within 
each subdomain, the data is represented pseudospectrally 
through an expansion into Chebyshev polynomials in the 
radial direction and spherical harmonics in the angular 
directions. When we take spatial derivatives, we simply 
use the analytical expressions for the derivatives of the 
individual basis functions and sum the results up using 
the expansion coefficients. The code is pseudospectral 
and not fully spectral in that we do not evolve the series 
expansion coefficients, but instead keep the values of E 
and B on a set of collocation points (at the intersections 
of the black lines in Fig.[^ optimized for translating back 
and forth into expansion coefficients (so that we can go 
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quickly into the series expansion representation and take 
spatial derivatives there, before jumping back). This way, 
we can implement the evolution equations in their natu¬ 
ral spacetime form, yet still take advantage of the high 
accuracy of spectral derivatives. The number of colloca¬ 
tion points represents the resolution of the simulation and 
corresponds to the highest order of the basis functions 
used (e.g. the largest I in We label the different 

levels of resolution using an integer k which changes lin¬ 
early with the total number of collocation points in each 
spatial dimension. In the radial direction, the number 
of collocation points is given hy k + 6. In the spherical 
directions, the highest I in the harmonics used is 2k + 7, 
which translates into 2fc -I- 8 collocation points in the 9 
direction, and Ak + 16 collocation points in the (j) direc¬ 
tion. To give readers an intuitive feel of the density of 
collocation points. Fig. plots the grid for fc = 6. 


B. Initial data solver 

It is frequently the case with FFE evolutions that the 
initial B held is not divergence free, and subsequently 
V B is cleaned using some additional cleaning held 
(see Appendix during evolution. We do not imple¬ 
ment such a cleaning held, and instead properly solve 
the constraints for our initial data. This is only neces¬ 
sary for the evolution of the perturbed solutions carried 
out in Sec. EYi and is not used for the numerical tests of 
Sec, [me] where we simply use exact constrain-satisfying 
analytical solutions as initial data. 

The divergence V • B can be removed by solving the 
Poisson equation 

= -V • B (35) 



Time / R 

+ 

FIG. 4: Top: The L 2 norm of the error measure for 

the constant B held simulation. Bottom: The L 2 norm of 
the error measure AE^. The simulation time is shown in 
units of the light-crossing time R+ from the origin to the 
outer boundary R+. There is only one curve being plotted in 
each panel, which oscillates quickly giving the appearance of 
a band. 


on the initial spatial hypersurface, and then B -b V<I> will 
be a divergence-free field. We solve Eq. (351 with the 
multidomain spectral method described in Ref. [89], and 
set the Dirichlet boundary condition $ = 0, so as to 
preserve the original B as much as possible by avoiding 
any unnecessary V $ pointing between different segments 
of the boundaries. 


We also tune the E field so that the FFE constraint 
E B = 0 is also satisfied by the initial data. This is 
achieved easily through an algebraic operation 


E -)► E - 


E B 

R2 


B. 


(36) 


As this only modifies the E field, it will not interfere 
with the earlier divergence cleaning step. We emphasize 
that this is the same operation as is given by Eq. (15) of 
Ref. [15] , except we are strictly applying the operation on 
the initial data, whereas in m this operation is applied 
at every time step of the evolution. 


C. Numerical tests 


1. Constant B field in flat spacetime 


A particularly simple solution to the FEE equations is 
a constant B field with vanishing E field in a Minkowski 
background spacetime. We choose the spatial slicings 
such that K = 0, and adopt a Cartesian coordinate 
system in which N = 1, f3 = 0, and g is the unit 
matrix. This is a physically trivial yet numerically in¬ 
teresting solution. Here and for the rest of the paper, 
the computational domain is a spherical shell, which is 
broken down into “subdomains” of concentric spherical 
shells (four shells extending from a radius of i?_ = 1.9 
to i?+ = 15 code units for this test; the speed of light 
is unity in the code units). Therefore, the constant B 
field subtends all possible angles with the normal n of 
the subdomain boundaries, including the special case of 
(n • B)^ = (n X E)2 = 0, when the evolution system 
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prescribed by Eqs. (15) and (161 becomes ill posed. 

For initial conditions, we use an analytical constant B 
field of unit strength without any added noise or per¬ 
turbation, which is also imposed as a Dirichlet bound¬ 
ary condition on the incoming (into the computational 
domain) characteristic modes during the evolution, on 
the external boundaries of the computational domain (at 
i?_ = 1.9 and at R+ = 15). The interior of the com¬ 
putational domain is evolved with Eqs. (15) and (16), 
and the differences between the numerical and analytical 
values of B and E provide a precise measurement of the 
simulation error. Letting AB and AE be the differences 
between the numerical and analytical values of B and E, 
we compute the L 2 norms of AB • AB and AE • AE: 


L2{s) 



dV = \/|det(h)|da;^ , 


(37) 


where h is the spatial metric and E is the computational 
domain. The values of L 2 {AB^) and L 2 {AE^) over the 
entire computational domain are plotted in Fig. From 
this figure, we see that despite not being strictly strongly 
hyperbolic, the evolution system can be evolved stably 
for a long period of time. 

The constant B field test is, however, not suitable for 
examining the convergence behavior of our pseudospec- 
tral code with increasing resolution, because the spatial 
derivatives of the evolution variables vanish identically 
in approximations to the spatial derivatives at any or¬ 
der. For this task, we turn to the nontrivial analytical 
solutions given in Ref. m- Once again, we utilize ana¬ 
lytical solutions because they provide precise references 
for comparison, allowing for more rigorous convergence 
tests. 


2. Analytical nulE solutions 

For a nontrivial analytical solution that’s more suitable 
for testing convergence, we turn to the time-dependent 
fully three-dimensional wave given at the end of Sec. |II Cj 
whose structure is plotted in Fig. This null'*' wave 
solution has a large amplitude and is fully nonlinear. 

In Fig. we plot the L 2 norms of the errors in B 
and E for our simulations of the null+ solution. Due 
to constraints on computational resources, we evolve the 
simulations to lOOOM, which is around 16 cycles for the 
time-dependent F{iy). We carry out evolutions at ten 
different resolutions where the number of radial colloca¬ 
tion points is given by fc -I- 6, fc € {0, • • • , 9}. Note our 
maximum Imax for the higher resolutions becomes exces¬ 
sive and the improvements in accuracy come mainly from 
the increase in radial resolution, but we keep the same 
Imax vs. k relationship to ensure consistency across all 
resolutions. Recall from Sec. |Tng the number of col¬ 
location points in each spatial dimension scales linearly 
with k. On the other hand. Fig. shows that the errors 
decline approximately exponentially with k, so our FFE 
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FIG. 5: Top: The L 2 norm of the error measure AB^ for 
the null’*' simulations. Bottom: The L 2 norm of the error 
measure AE^. There are ten resolutions being plotted, with 
the number of radial collocation points given by A: -|- 6, and 
the maximum I for the angular y*"* decomposition given by 
2k + 7. In both panels, the error decays with increasing k, 
so that k = 0 corresponds to the topmost (black) line, and 
the lower lines are, in turn, k — 1,2,3, ■■■ . The errors decline 
approximately linearly in log scale, indicating that our pseu- 
dospectral code achieves exponential convergence as expected 
over a wide range of resolutions, with the highest resolution 
being limited in accuracy by machine precision. 



implementation achieves exponential convergence for this 
time-dependent, nontrivial, but smooth (without current 
sheets) test case, as expected from its pseudospectral na¬ 
ture. 


3. Quasinormal modes of the black hole magnetosphere 

Next, we present a test case with richer dynamics. 
The test is designed to demonstrate the code’s ability 
to correctly treat various types of waves allowed by the 
force-free equations, and in particular account for their 
interactions with the spacetime curvature to a high ac¬ 
curacy. As mentioned in the Introduction section, the 
null+ solutions we examine are envisaged to play the role 
of carrying energy across magnetospheres of black holes 
or neutron stars, therefore a more direct demonstration 
of the relativistic effects is beneficial. 
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Specifically, we launch a train of ingoing FFE waves 
towards the black hole, and examine the magnetosphere 
quasinormal modes (QNMs) being excited. Similar 
modes are intensely studied for gravitational waves, and 
constitute the main postmerger signal after a coales¬ 
cence of two black holes. More recently, they have been 
shown to also be present in a black hole magnetosphere, 
where the so-called trapped FFE modes (in the short- 
wavelength/eikonal limit, they become the fast magne- 
tosonic waves, so they will be referred to as generalized 
magnetosonic waves below) can be excited and trapped 
by the gravitational potential well of the black hole, and 
leak out over time with an exponentially decaying ampli¬ 
tude. The frequency and decay rates of these modes are 
completely determined by the spacetime geometry and 
the nature of the wave, and is independent of how they 
are excited. In particular, this means whatever the fre¬ 
quency of the wave train we send in, we should observe 
the same frequency and decay rate for the QNM it ex¬ 
cites. 


When we launch a wave train consisting of in general 
the long wavelength generalizations of both Alfven and 
magnetosonic waves towards the black hole, part of the 
wave is swallowed by the black hole, while the rest scat¬ 
ters around it and travels back out. During the pro¬ 
cess, some generalized magnetosonic waves are excited 
as QNMs, which do not have a counterpart consisting of 
generalized Alfven waves. Our code needs to be able to 
correctly distinguish the propagation behaviour of both 
types of waves, especially their different interactions with 
spacetime curvature, in order to reproduce the frequency 
and decay rates predicted by analytical calculations. 


The wave train is constructed by first building a 4>q 
(the NP scalar representing ingoing waves) that has a 
spherical harmonic angular profile and a sinusoidal radial 
profile Asin(D(r—rg)) with rg = 200M, enclosed in a top- 
hat-like radial envelope extending from lOOM to 300M. 
We also add & q = 1000 magnetic monopole which enters 
through the Coulomb background as in Eq. ( [^ . To¬ 
gether, these scalars allow us the build the Faraday tensor 
via Eq. ( [33| ) and then E and B via Eq. (34) (Fig. [^pro¬ 
vides a visual depiction of the wave train and its time evo¬ 
lution) . In order to accommodate the wave train, we use 


a slightly different domain structure from that of Sec. IIC 


and the rest of the paper (including the stability studies 
later). Namely we have 16 concentric spherical shells ex¬ 
tending from 1.05M to 495M, with the number of radial 
collocation points given by 5 -I- 2K and those for the 0 
and (j) directions by 7 -I- and 14 -|- 2K (we have used 
capital K to dist inguish it from our usual resolution label 
of k used in Sec. 111C2 and later stability test sections). 

The bottom panel in Fig. shows both the power and 
the limitation of the spectral methods. The gray surface 
in that figure is constructed by linearly interpolating be¬ 
tween the collocation points (done by the visualization 
software, and is not the actual spectral representation 
of 3?<(>g), which is not smooth at all because there are 
very few collocation points inside each radial oscillation 


Time: 0 


Time: 100 



FIG. 6: Top row: Contours of that show the structure of 
the wave train. It has a angular profile, and consists of 
fast radial oscillations with frequency Q. It is confined into 
a radial top-hat-style envelope of width 200M. The shape of 
the wave train shows that our test is fully three dimensional, 
without any axisymmetry. The figures on the left and right 
depict the wave train at earlier and later times, showing that 
it is initially moving inwards toward the black hole. The green 
semitransparent surface is the equatorial slice of the computa¬ 
tional domain, included to indicate its extent. The grid struc¬ 
ture on this slice is also plotted to provide a visual description 
of the density of collocation points at A = 7. Bottom figure: 
The ?R.(fio is depicted by warping half of the equatorial slice 
with a vertical displacement proportional to the ^(j>o values 
at collocation points, and then connecting them via “straight 
lines” into the gray surface (i.e. the gray surface is a trivial 
interpolation between the collocation points and is not the 
actual 5R(/)o as represented by the spectral functions). The 
green surface is the tophat-like radial envelope function. This 
figure also corresponds to A = 7. 


period. The actual spectral “interpolation” using our 
Chebyshev polynomials on the other hand, gives very 
smooth functions (see Figs. [^[^ andbelow). This shows 
that the spectral methods are capable of providing accu¬ 
rate representations of functions with very few colloca¬ 
tion grid points, in other word their effective resolution, 
given a fixed number of grid points, is very high. On the 
other hand, the black square in that figure highlights the 
region close to the rather sharp corner of the top-hat-like 
envelope of the form 



(38) 


with r_|_ and r_ being the outer and inner boundaries of 
the envelope respectively. The spectral methods behave 
rather differently from finite difference methods, where 
the sharp kinks are simply under-resolved and rounded 
due to numerical dissipation. With the present low radial 
resolution, the steep envelope cannot be resolved, and 
leads to Gibbs oscillations bleeding into the surrounding 
region. While the shape of the waveform is thus not per¬ 
fectly represented, we note that nevertheless the shape 
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FIG. 7: The log of 4'^’^'^ plotted as a time series. The curves 
corresponding to different wave train frequencies If are shifted 
slightly relative to each other vertically, so that the curves 
won’t lie right on top of each other, and we can see more 
clearly. The flat region for the curves correspond to the 
wave train itself, moving back outwards after having scat¬ 
tered around the black hole. The sloped region corresponds 
to the exponentially decaying QNM. We can see that the wave 
train has very different frequencies, but the QNMs excited in 
all cases share the same frequency and decay rate. 



Analytical Ref. 

n = 0.2 

n = 0.3 

II 

0 

n = 0.5 

UJ 

0.457596 

0.444071 

0.453783 

0.453742 

0.44466 

7 

0.0950044 

0.0922016 

0.0885479 

0.0893736 

0.0901566 


TABLE I: The fitted values of the frequency lu and decay rate 
7 of the QNMs excited by wave trains of different fl, as well 
as the reference analytically computed values. 


is propagated without diffusive broadening [5D] (also see 
this reference for an additional suite of tests on the nu¬ 
merical behaviour of the pseudospectral infrastructure 
underlying our force-free code), as typically caused by 
numerical viscosity of finite-difference codes. 

To examine the QNMs, we need to extract the out¬ 
going wave component (j >2 of the Faraday tensor (QNMs 
escaping from the trapping gravitational potential will 
show up as an exponential tail to the wave train in ^ 2 )- 
We do so by extracting the (j )2 values via Eq. (21) on a 
sphere located at lOOM from the coordinate origin, and 
integrate this <p2{9, (f>) distribution against the spher¬ 
ical harmonic to get a single scalar value representing the 
(2, 2) harmonic component of (f> 2 - We do so at each time 
step and obtain a time series for this harmonic coefficient 
of the outgoing wave, which we denote ^2 ’ We note 
that in order to make comparison with perturbation the¬ 
ory that predicts the QNM frequencies, we use a small 
wave amplitude oi A = 0.0001. There is however still 
nonlinear interactions between the wave and the back¬ 
ground magnetic monopole. 

We first test whether the correct QNM wave is pro¬ 
duced in our code, specifically, whether it has the correct 
frequency and decay rate as predicted by analytical mod¬ 
els. To this end, recall that these values are independent 



FIG. 8: The dashed black lines indicate the htting results 
(they extend only over the fitting interval) corresponding to 
the fitted parameters in Table |I] while the curves being htted 
are the negatively sloped segments of the curves in Fig. 


of the specifics of the initial wave train, so a stringent 
test consists of launching initial wave trains with differ¬ 
ent frequencies Q, and see if we get the same set of QNM 
parameters despite this major difference. If we do obtain 
the same parameters as expected, we can then be confi¬ 
dent that the code has managed to cleanly excite and ex¬ 
tract the QNMs, and the decaying tail is not some other 
leftover feature from the wave train (such as the oscilla¬ 
tions seen in the black square of Fig. which is in reality 
much smaller in amplitude and drops off faster than the 
QNMs). We simulate four cases with Q = 0.2,0.3,0.4, 
and 0.5 at K = 7, and plot their in a log plot 

(so exponential decay shows up as a straight line with 
a negative slope equaling the decay rate) in Fig. 0 By 
visual inspection, we can already see that despite the 
rather different frequencies of the wave train, the fre¬ 
quencies u! and decay rates 7 of the exponential tails in 
the four cases are indeed the same. We can also carry 
out a more quantitative fitting of these quantities using 
Mathematica’s FindFit function. The results and the an¬ 
alytical reference values are shown in Table. |T] (We also 
plot in Fig.|^the fitting results (dashed black lines) inside 
the fitting interval to provide a visual depiction of the fit¬ 
ting quality. ) We indeed find a good match between the 
measured and predicted values. Aside from truncation 
error, the excitation of QNMs of other overtone numbers 
n > 1 (we have only considered the fundamental overtone 
n = 1 that decays the slowest) is likely a main source of 
any residual errors. 

The second test we carry out is a convergence test, 
to verify that our code still achieves exponential conver¬ 
gence with this dynamically richer setup. To this end, 
we fix Q = 0.4 and perform four simulations at K = 5, 
K = 7, K = 9 and K = 11. The log plots of for 
these runs and their differences are shown in Fig.[^ which 
shows that exponential convergence is indeed achieved. 
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FIG. 9: The top figure shows the wave train and the QNM tail 
for the four resolutions of the Q — 0.4 simulation. The bottom 
figure shows the differences between resolutions. Because the 
differences are approximately equally displaced in the vertical 
direction, which is plotted in log scale, the simulation indeed 
achieves exponential convergence. 


IV. STABILITY OF THE NULL+ SOLUTIONS 


In this section, we utilize our new FFE code and carry 
out numerical simulations in order to determine whether 
the null'*' solutions are stable, in the sense of whether a 
perturbed solution will asymptote to an exact null’^ so¬ 
lution over time, or diverge from it. As an underlying 
unperturbed null'*' solution, we use the same configura¬ 
tion shown in Sec. m and studied in Sec. C 2[ with 
q — 1000, F{l') = Acos{fli') and A = 1. However, here, 
we will vary H. Note that we have chosen a large absolute 
magnitude for the Faraday tensor in order to carry out 
a more stringent test, as small magnitudes will dimin¬ 
ish the significance of the nonlinear terms. We also use 
the same domain structure as in Secs. |II C and III C 2 


namely eight spherical shells extending from 1.95M to 
195M. 


We denote the Faraday tensor for the unperturbed so¬ 
lution as F and impose it as a Dirichlet boundary condi¬ 
tions on the incoming characteristic modes as discussed 
in Sec. |III A| In other words, we enforce the condition 
that there are no incoming perturbative modes. On the 
other hand, there is no restriction on the outgoing char¬ 


acteristic modes, so overall, we have a purely outgoing 
boundary condition for the perturbations, just as when 
one studies the mode stability of the Kerr metric by cal¬ 
culating its QNMs, for example. 


A. Perturbation in (f>o 


The simplest perturbation to the null+ solution is a 
variation in the initial (j)o profile, so that it no longer sat¬ 
isfies Eq. (29). Specifically, we generate a Faraday field 
from an altered (j)Q (in addition to an unperturbed 
monopole component) that is perturbed away from the 
exact solution (Eq. 29) by 


, sm 


S4>o = 0.25—— F{v)i sin'f/j, 

Ap 


(39) 


where we choose the same Flv) function as the unper¬ 
turbed solution. For simplicity, we also choose = 0 so 
that the unperturbed background solution is time inde¬ 
pendent. Such a solution does not have the usual char¬ 
acter that one would associate with a travelling wave, 
although it is still technically a wave in the same sense 
that a constant would satisfy a simple 1-D wave equa¬ 
tion. Nonetheless, this solution is still physically inter¬ 
esting, as the energy flux does not vanish, so that the 
solution describes an “electromagnetic wind” or “Poynt- 
ing wind” m- We will examine a time-dependent case 
later in Sec. II V Cl 

Given the altered field c/jq, we now define our perturbed 
initial data. Because we employ boundary conditions 
based on the unperturbed background solution, we re¬ 
quire the perturbed initial data F^ to approach the un¬ 
perturbed solution at the boundaries (strictly speaking 
we don’t need a boundary condition at the inner bound¬ 
ary as it is inside the event horizon, but we use the Dirich¬ 
let condition there for simplicity). We achieve this by 




FIG. 10: Left: The right-hand side of Eq. (35), or explicitly, 
the initial — V-B value on a vertical slice of the computational 
domain before we apply our initial data solver (the z-axis 
corresponds to 0 = 0 and 6 — n). Right: The left-hand side 
of Eq. (35), which is the value of that emerges after 
solving for the initial data. 
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FIG. 11: Top: The I /2 norm of the difference measure 
for the nulh^ simulations initially perturbed in the i^o values. 
Bottom: The L 2 norm of the difference measure AE^. The 
arrangement of the collocation points is the same as for Fig.[^ 


blending between and F via 

F^ = /F^ + (1-/)F, 
where the weighting function is 

f = P^(2-PfQ^(2-Qf 


with 


P = 2 -^—and Q = 2- . 
R+ — R- TT 


(40) 

(41) 

(42) 


Equation (41) ensures 0 < / < 1 such that the unper¬ 


turbed solution dominates on the domain boundaries, 
ensuring a smooth transition to the Dirichlet boundary 
conditions imposed there. The angular dependence in A 
ensures that the unperturbed solution also dominates on 
the vertical axis (0 = 0 and 9 = tt). This feature is not 
constructed to satisfy any requirements of the present 
form of perturbation, but instead is included for later 
convenience in Sec. HYH 

The perturbed solution does not automatically satisfy 
the constraints; since the monopole resides inside the 
event horizon and outside of our computational domain, 
the magnetic field should be divergence-free within the 
computational domain. We are therefore free to solve 


the Poisson equation in Sec. |IIIB| to remove any diver¬ 
gence. In Fig. 10 we plot the right and left-hand sides of 
Eq. (351 and show that the initial data solver performs 
as designed by removing V • B. Note that there is some 
high-frequency noise in the output of this elliptic solver 
(see the center of Fig. 10 (b)), which we partially remove 
by passing the initial data through a spectral filter that 
reduces the high-frequency spectral coefficients in the ra¬ 
dial direction. The removal is only partial as the filtering 
strength is chosen to be conservative so it does not intro¬ 
duce new divergence into the magnetic field. As a second 
step, we also use Eq. (36) to impose the FEE constraint 
E B = 0. 

we plot the evolution of L 2 {AB^) and 


11 


In Fig. 

L2{AE^), ^ere AB = B^-B is the difference between 
the evolved B^ and its unperturbed counterpart B (com¬ 
puted analytically). AE is defined similarly. We observe 
that the differences drop by several orders of magnitudes 



FIG. 12: Top four panels: The relative difference | AB|/|B| on 
a vertical slice of the computational domain. The magnitude 
of the relative difference is indicated by the height of the sur¬ 
face as well as the colouring. Different panels correspond to 
different times. The initial perturbation seen at t = 0 prop¬ 
agates inwards creating the pattern shown at f = 60. After 
scattering around the black hole, the perturbations propagate 
outwards, forming the pattern seen at t = 150, and eventually 
begin to exit the computational domain as seen at t = 270. 
The bottom two panels show the initial and later perturba¬ 
tions on the horizontal equatorial slice of the computational 
domain, they show that the perturbations are three dimen¬ 
sional in nature, without axisymmetry. 
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FIG. 13: The numerical error in B** (as seen in |ABl/lBi) 
being generated at subdomain boundaries (signified by dense 
concentrations of black lines). The gray frames highlight the 
locations of the error on such boundaries. Note that the warp¬ 
ing scale and the colour map are different from Fig. |12[ b ut 
the plane shown is the same as those appearing in Fig. |12| 


as time progresses. A more detailed distribution of AB 
is shown in Fig. which plots the relative difference 
|AB|/|B| on a vertical slice of the computational domain. 
The magnitude of the relative difference is indicated by 
the height of the surface as well as the colouring. The 
panels correspond to different times and suggest that the 
perturbation does not diverge, but instead propagates in¬ 
wards to begin with, and then outwards after scattering 
around the black hole, before exiting through the outer 
boundary. 


Consequently, our observation indicates that there are 
no diverging modes excited by our initial perturbation 
at the fully nonlinear level, and so the null’*' solution is 
stable against our perturbation. Furthermore, the even¬ 
tual exit of the perturbation is consistent with the null'*' 
solution being asymptotically stable (attracting). We 
note however that we cannot rigorously prove the lat¬ 
ter stronger stability, as L2(Ai3^) and L 2 {AE‘^) values 
at late times do not reach their small sizes in the un¬ 
perturbed convergence test case shown in Fig. The 
main culprit appears to be the high-frequency noise in 
the initial data we see in Fig. 10 (b), which is absent from 
the analytical initial data used for the convergence tests. 
The fact that the convergence behaviour improves after 
we apply spectral filtering to the initial data before start¬ 
ing the evolutions provides evidence for this conclusion. 
In particular, the SpEC code utilizes a penalty method 
[HSHIHI to enforce consistency across internal boundaries, 
which allows for discontinuities to exist temporarily. It 
has been noted in several previous studies [Hniini] that 
high-frequency noise tends to induce large discontinuities 
at the internal boundaries, thereby creating errors in B. 
The situation is the same for our FFE evolution system, 
as can be seen in Fig. The high-frequency noise also 
destroys the convergence at the highest resolutions (see 


X 

FIG. 14: Example congruence of geodesics in the equatorial 
plane of the Schwarzschild black hole, with a nonvanishing 
impact parameter b translating into V = 2/27M^ (not the 
value we use for the simulations, but instead chosen to ac¬ 
centuate the perturbation). The red circle is the location of 
the event horizon. Angle Y is the angle at which the dashed 
geodesic strikes the origin. 


Fig. 11), because lower resolution (smaller k) acts as an 
effective spectral filter, shielding the lower resolution sim¬ 
ulations from high-frequency noise. We expect this com¬ 
plication to disappear as we ascertain the source of the 
high-frequency noise (one possibility is the boundary con¬ 
dition of d* = 0 being too simplistic) and further improve 
our procedure of solving Eq. (35). 


B. Perturbation in the propagation direction 

With numerical experiments, we can only state that 
there are no diverging modes being excited by the par¬ 
ticular perturbations that we introduce into the initial 
data. Therefore, in order to provide as strong evidence 
as possible for stability (meaning the nonexistence of di¬ 
verging modes in general), it is important that we con¬ 
sider a set of initial perturbations that is as general as 
possible. In the last section, although we started with 
specific modifications to (/jq, our initial data construction 
procedure nevertheless has to go through the blending 
and constraint solving stages, which means the resulting 
perturbation is in fact rather general (containing many 
modes). 

However, the wave component (the (jjQ piece of the 
Earaday tensor) of the simple perturbation studied in 
Sec. |IV A| still follows the GPNDs initially as we have 
retained the use of the Kinnersley tetrad when construct¬ 
ing the perturbed Faraday tensor (^o is by definition the 
wave propagating along the tetrad’s ingoing null direc¬ 
tion, which for the Kinnersley tetrad is in a doubly degen¬ 
erate GPND direction) and is therefore (at least initially) 
not necessarily backscattered by the spacetime curvature 
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(see Sec. |V]). In other words, the perturbation may be 
restricted to a rather specific submanifold in the solu¬ 
tion space. In principle, it is possible that even though 
null"'" solutions are stable on this submanifold, they may 
still be unstable under perturbations off of it. Although 
we would expect numerical truncation errors to induce 
perturbations away from this submanifold regardless of 
whether the initial perturbation places us on it or not, the 
time scale for this occurring may simply be too large so 
we fail to observe the possible instabilities in the previous 
test, and a more appropriate off-the-submanifold initial 
perturbation is desirable. In this subsection, we con¬ 
struct perturbations that are nontrivial in the sense that 
the wave component of the initial data does not follow de¬ 
generate GPNDs everywhere and is therefore generically 
backscattered by spacetime curvature. In other words, 
the initial data do not possess one of the core features of 
the null"'' solutions, and the perturbation is not restricted 
to some specific submanifold in the solution space. 

For the initial wave propagation direction, we use 
a congruence of null geodesics (that is generally not 
tangential to degenerate GPNDs) in the Schwarzschild 
spacetime, which admits analytical descriptions. We 
will begin by introducing the procedure for generat¬ 
ing a single geodesic within the congruence, building 
an adapted Newman-Penrose null tetrad along it whose 
ingoing null basis vector is tangential to the geodesic, 
and constructing the field whose wave component 
4>q (as defined with respect to that newly built tetrad) 
follows the geodesic direction. Later, we will describe 
how to obtain the entire congruence, thus filling in the 
field everywhere. We begin by choosing Boyer- 
Lindquist/Schwarzschild coordinates {r,d,(j)) so that the 
geodesic lies on the equatorial plane. The null geodesic 
then satisfies the equation 

O "V-929-93, (43) 

where 

Ml 1 1 ( mV ^ 



FIG. 15: Similar to Fig. |10[ but for initial data with a per¬ 
turbed propagation direction. Left: The right-hand side of 
Eq. ( |35[ ). Right: The left-hand side of Eq. ( |35[ ). 



Time / M 


FIG. 16: Top: The L 2 norm of the difference measure for 
the null’*’ simulations initially perturbed in the propagation 
direction. Bottom: The L 2 norm of the difference measure 
AE^. 


with V — 1/b^ for impact parameter b. With small im¬ 
pact parameter {V > 1/27M^), the null geodesic will be 
absorbed by the black hole, and the shear-free principal 
congruence tangential to the degenerate pair of incoming 
GPNDs corresponds to 6 = 0. The solution to Eq. (43) 
is then given by y(4>) = p{(j) -f F|( 72 , 5 ' 3 ), where p is the 
Weierstrass elliptic function and ( 52 , 53 ) are its invari¬ 
ants. The angle Y is the angle at which the geodesic 
strikes the origin r = 0. 


Armed with the geodesic, we can now build a Newman- 
Penrose null tetrad adapted to it. First, we calculate the 
spatial tangent to the geodesic, and then we convert it 
into a null 4-vector and apply the Jacobian from Boyer- 
Lindquist to Kerr-Schild coordinates (see Appendix[D|for 
details). We will let the n basis vector be in the direction 
of this four-dimensional null tangent, while keeping 1 the 
same as that of the Kinnersley tetrad. The scaling of n 
is then fixed by Ena = —1- The remaining m and m 
bases can be fixed using a Gram-Schmidt procedure [S^ . 
We first define two spatial vectors C and D, so that m = 
l/y2(C-biD). Let G and H be the Kinnersley tetrad 
versions of C and D, respectively. We can then achieve 
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proper orthonormality for the new tetrad by setting 


G“ = G“ + G'’4n“ + G''’nfar, G“ = 


G“ 


VG^Gi 


, (45) 


and then 


- WGbG'^, (46) 

fja 


£)“ = 


\/HWb' 


(47) 


Under this tetrad, we choose a preliminary Faraday field 
F'^ according to Eq. (22), with the (j)o distribution pre¬ 
scribed by Eq. (29), which now describes an incoming 


wave (there is no ((> 2 ) travelling in the new n direction 
that is different from the GPNDs provided b ^ 0. 

We now turn to filling the entire computational do¬ 
main with geodesics and the Faraday field tensor. To 
this end, we begin by specifying an impact parameter b 
and populate the equatorial plane with null geodesics by 
varying Y (see Fig. 14). We then take the a; > 0 portion 


of Fig.pT|and rotate it around the z-axis, thus filling the 
entire 3D space. However, when the impact parameter 
does not vanish (for our perturbative study, we choose 
b — M/i/Id), the resulting congruence will be singular 
on the vertical axis. We eliminate this problem by con¬ 
structing an unperturbed null solution with the same 4>q 
(but with & = 0), and blend its Faraday tensor F with 
the F'^ associated with the b ^ 0 congruence, so that 
only the unperturbed null solution F that is regular on 
the vertical axis is present there. Explicitly, we set 


F^ = /F^ + (1-/)F, 


(48) 


Time: 0 


Time: 60 





FIG. 17: The relative difference | AB|/|B| on a vertical slice of 
the computational domain. The initial perturbation seen at 
t = 0 propagates mostly outwards, creating the patterns seen 
at t = 60 and t = 100, and exits the computational domain 
as seen at t = 200. 



FIG. 18: Similar to Fig. |15| but for initial data with a per¬ 
turbed propagation direction based on a time-dependent un¬ 
perturbed solution. Left: The right-hand side of Eq. (351 
Right: The left-hand side of Eq. (|35 


with the weighting function / given by Eqs. (41) and 
(42) (the choice of Q in those equations was made 


in anticipation of this regularization procedure on the 
poles). Furthermore, we will add an unperturbed mag¬ 
netic monopole with q = 1000 t o obtain a perturbed 
null+ solution just as in Sec. IV A We also carry out the 
same constraint enforcement procedure as in Sec. |IV A| 
For this section we will keep the background solution time 
independent with D = 0, and leave the time-dependent 
case to the next section. 

We note that just as before, the blending and con¬ 
straint solving stages ensure that the perturbation we 
start the simulation with is in fact rather general, with 
many radial and spherical modes excited. Our procedure 
differs from an explicit sum of modes with random co¬ 
efficients in that the specific modifications to the wave 
propagation directions essentially introduce correlations 
into the mode coefficients, so that the modes don’t acci¬ 
dentally cancel out, leaving us with a perturbation with¬ 
out propagation direction change. It also helps to avoid 
the initial data solver simply removing the uncorrelated 
constraint-violating modes and bringing us back to the 
unperturbed exact solution. 

The output of the initial data solver is displayed in 
Fig. while the evolution of AB and AE is shown 
in Figs. [ 1 ^ and Despite having a different type of 


initial perturbation, we observe a similar behaviour as 
seen in Sec. |IV A[ with no diverging modes occurring, 
and with the perturbation eventually exiting through the 
outer boundary. 


C. time-dependent background solution 

For our final numerical setup, we introduce a time de¬ 
pendence with n — 0.1, so that the background solution 
has the familiar character of a travelling wave. The pro¬ 
cedure for introducing perturbations is otherwise identi¬ 
cal to Sec. |IVB[ so that the waves are initially travelling 
in directions different from the degenerate GPNDs. 

The output of the initial data solver and the evolu- 

Despite 


tion codes are displayed in Figs. 18 19 and 20 
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the change of energy flux character from electromagnetic 
winds to waves, our simulation suggests that the now 
time-dependent null+ solutions are also stable. The most 
noticeable difference with the two earlier cases is that the 
perturbation propagates almost entirely inwards initially, 
and is absorbed by the black hole almost completely after 
one light-crossing time. 


V. SCATTERING BY SPACETIME 
CURVATURE AND THE ROLE OF GPNDS 

In this section, we seek to shed some light on the ques¬ 
tion of what feature of the analytical solutions examined 
in this work allows them to avoid being backscattered 
by spacetime curvature. Because we will only carry out 
analytical studies on the unperturbed exact solutions in 
this section, the null and null+ solutions are exactly the 
same in terms of their wave propagation properties. So 
we will consider pure null solutions for brevity, with the 
understanding that the conclusions translate to the wave 
component of the null^ solutions trivially. 

The answer to the scatter-avoidance question is inter¬ 
esting in that it provides guidance to the search for sim¬ 
ilar FFE solutions in other spacetime backgrounds, or 
solutions to non-FFE equations. For example, when sev- 



Time: 0 


Time: 70 



Time: 200 



FIG. 20: The absolute difference |AB| on a vertical slice of 
the computational domain. The initial perturbation seen at 
t = 0 propagates inwards creating the patterns seen at t = 
70 and t = 120. By t = 200, the perturbation has been 
almost entirely absorbed by the black hole. We have shown 
the absolute, rather than the relative difference, because |B| 
increases quickly when we approach the black hole, so it is 
more difficult to see what is going on in the inner regions 
with a relative difference plot. 
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FIG. 19: Top: The I /2 norm of the difference measure AB^ 
for the initially perturbed time-dependent null'*’ simulation. 
Bottom: The L 2 norm of the difference measure AE^. 


eral analytical solutions mm to the FFE equations were 
first found, it was not immediately clear why such simple 
solutions exist nnnn], given that the EFE equations are 
nonlinear. Furthermore, such scatterless null solutions 
are closely related to important advances in mathemati¬ 
cal physics, such as the discovery of new solutions to the 
Einstein equations [53 ES] , and the definition of twistors 
m- Therefore it is informative to try to understand the 
core features of these solutions at an intuitive level. 

A hint on the answer to this question is provided by the 
Goldberg-Sachs theorem [98], which states that scatter¬ 
avoiding waves must propagate along a repeated princi¬ 
pal null direction (GPND"*) of the Weyl curvature tensor. 
However, as far as the authors are aware, there is no ex¬ 
plicit analysis in the literature of the reverse question, 
i.e. whether all waves propagating along GPNDs are to 
some extent scatter avoiding. Aside from shedding some 
light on the scatter avoidance puzzle, this analysis also 
fills a gap in the literature by providing a simple physical 
intuition on the concept of GPNDs, which underlies such 
important constructs and results as the Petrov classifica¬ 
tion of spacetimes and the peeling theorem [^[99HT02]. 

Before diving into the technical details, we first sum- 


We denote these principal null directions as GPNDs, rather than 
PNDs, to emphasize that they are gravitational PNDs, and dis¬ 
tinguish them from the electromagnetic PNDs that we will en¬ 
counter later. 
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marize the results of this section, and provide some intu¬ 
itions as to why one should expect these results. Readers 
only interested in the conclusions and their uses can skip 
the derivations presented later in this section. The con¬ 
clusions are 

1. Null solutions that propagate along repeated GP- 
NDs of multiplicity 3 or above will not experience 
scattering by spacetime curvature at all. This re¬ 
quires the background spacetime to be of Petrov 
type III or type N [10311104j (and of course type O 
which is flat) and for the null waves to be prop¬ 
agating along the special direction in those space- 
times identified by the repeated GPNDs. 

Such a situation arises when analysing the coinci¬ 
dent gravitational and electromagnetic wave signals 
from distant sources in the context of multimessen¬ 
ger astronomy. The two types of waves travel in 
the same direction for the majority of their jour¬ 
ney, and the gravitational wave gives a metric per¬ 
turbation of type N with the four-fold degenerate 
GPND pointing in the propagation direction. This 
means one does not need to worry about the grav¬ 
itational wave changing the electromagnetic wave 
during their long journey to Earth. 

2. When the spacetime is less special in the sense 
that the repeated GPNDs have less multiplicity, the 
null electromagnetic solutions travelling along the 
repeated GPNDs will in general experience some 
scattering by the spacetime curvature. However, 
the more severe backscattering (where a wave prop¬ 
agating in the opposite direction and/or a Goulomb 
background piece is created by the scattering) can 
be avoided, and the scattering only manifests itself 
as an influence on the wavefront profile. 

As Schwarzschild and Kerr spacetimes are of Petrov 
type D (the multiplicity of the degenerate GPNDs 
is two), this is the situation relevant for the exact 
analytical solutions (more precisely the null part of 
the null"*" solutions) presently under examination. 
The avoidance of backscattering allows for, e.g., 
setting (j )2 (the outgoing wave component of the 
Faraday tensor) and (j)i (the Coulomb background 
piece) to zero, and only solving for c/o (the ingo¬ 
ing wave piece) in the force-free equations, which 
significantly reduces the complexity of the solution 
finding process. We expect similar features to also 
be available when solving the equations of other 
field theories in these important spacetimes. 


® Petrov type I: four different GPNDs; type II: two degenerate 
GPNDs and two nondegenerate GPNDs; type D: two sets of 
doubly degenerate GPNDs; type III: a triply degenerate set of 
GPNDs and a nondegenerate one; type N: four-fold degeneracy 
in GPNDs; type O: Weyl curvature tensor vanishes. 




FIG. 21: The principal directions (black lines) and the asymp¬ 
totes (red lines) on curved surfaces. Figure on the left depicts 
the situation with a more generic surface, while the figure 
on the right depicts a more degenerate surface where two 
asymptotes become coincident and the principal direction in 
between them becomes flat. 


These conclusions create the obvious impression that 
the GPNDs are somehow particularly flat directions of 
spacetime, such that waves propagating along them see 
less of the spacetime curvature. This intuition can be 
made more visual by invoking an analogy with a curved 
surface embedded in a three dimensional (flat) ambient 
space. In fact, the name “principal null directions” al¬ 
ready invoke this analogy, but we will argue that instead 
of the principal directions, it is really the asymptotes on 
the curved surface that these GPNDs are more akin to. 
Given our quick scan of the literature, our discussion ap¬ 
pears to be the first to state this way of intuiting the 
GPNDs, and we hope it would be helpful to researchers 
newly acquainted with these important quantities. 


The curving shape of a surface embedded in the ambi¬ 
ent (see Fig. 21) is given by the extrinsic curvature tensor 
K, which is a rank two tensor. The twice contraction of 
a vector tangential to the surface at a point (the cen¬ 
ter/origin in Fig. 21) with K gives the curvature of a 
geodesic on the surface developed along that vector di¬ 
rection. The principal directions are defined to be those 
directions whose associated geodesics have the maximum 
or minimum (most negative/bending the other way) cur¬ 
vature among the different direction choices, and are in 
fact the eigenvectors of K. Their associated geodesics are 
the thick black lines on the surfaces in Fig. which are 
clearly the most curved (either in the up or down direc¬ 
tion) directions on the surface. The asymptotes on the 
other hand are defined to be those directions whose twice 
contraction with K vanishes. These are the red curves in 
Fig. 21 and are clearly locally flat/straight. Later in the 
section, we will show that the contraction of quantities 
representing GPNDs with those representing the space- 
time curvature tensor gives us zeros instead of maxima 
and minima, so the GPNDs are really more akin to the 
asymptotes. 


The principal directions are always half way in between 
the asymptotes, so these two sets of quantities are triv¬ 
ially related. When two asymptotes coincide, the princi¬ 
pal direction in between them is then forced to become 
coincident with both asymptotes and thus take up van¬ 
ishing curvature (the right panel of Fig. 21 shows this 
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situation, where there is also a thick black line aligned 
with the coinciding red lines at the bottom of the “val¬ 
ley”), so the spacetime becomes flatter in a sense as one 
of its extremally curved directions become flat (in the 
case of the curved surface, the surface on the right of 
Fig. that has coinciding asymptotes becomes devel¬ 
opable - its intrinsic curvature vanishes and the surface 
can be unfolded into a simple flat plane). The analog 
with this on the spacetime side is that more coinciding 
GPNDs signal that the spacetime is “flatter”, especially 
along the direction of those coinciding GPNDs (this is 
essentially the intuitive meaning of the Petrov classifica¬ 
tion). It is then not entirely surprising that FFE waves 
travelling along the “extra flat” degenerate GPNDs will 
experience less curvature-induced scattering. 

This analogy can be developed much further ® , but 
a detailed discussion will lead us too far on a digression 
away from the main content of this paper, so we stop 
here and turn to some derivations directly relevant for 
curvature-related scattering that back up the conclusions 
listed above, which also provide some more concrete ex¬ 
amples of the type of vanishing contractions underlying 
our analogy-based intuitive picture. 

The discussion below will rely heavily on the spinor for¬ 
malism, which reveals the characteristic structure of the 
Weyl curvature tensor in a significantly more transpar¬ 
ent manner than the tensor formalism. We include a brief 
summary of spinors in Appendix]^ The most important 
feature for us is that the spinors can roughly be seen 
as “square-roots” of null vectors with the tensor prod¬ 
uct of a pair of complex conjugate spinors 0"^ and 

0 A corresponding to a null vector. In addition, the self¬ 
contractions of the spinors vanish ( 0 a 0 ^ = 0 = 0 a' 0 ^ ) 
just as they do with null vectors. 

The null solutions of Ref. [I2] are the FFE counterparts 
to the vacuum null electromagnetic solutions described in 
Ref. m- A null solution is defined by the property that 
the two principal null directions (EPNDs^) of the Eara- 


° We brief mention why the name “principal null direction” is his¬ 
torically given to the GPNDs, even though they are more like 
asymptotes, and not the principal directions on a surface. That 
name assignment comes from their being related to some eigen¬ 
value problem of the spacetime curvature tensor when written 
in the tensor language. However, when we migrate to the spinor 
language, the eigenvalue problem switches into one over some 
other directions that sit in between the GPNDs (they are the 
null basis vectors of the so-called canonical transverse tetrads 
[571. and they are really the things that are akin to the princi¬ 
pal directions on a surface), just like the principal directions on 
a surface sit in between the asymptotes. So in the spinor lan¬ 
guage, the “true” nature of the GPNDs become more apparent, 
and using something called sectional curvatures 1105111061 , we 
can show that it is in the spinor and not the tensor language, 
that the eigenvalue problem is more closely related to the curva¬ 
tures of geodetic submanifolds of the spacetime (analogs to the 
curves on a surface), and so the spinor language is the one that’s 
more appropriate for building analogies. 

^ We will use EPND to refer to principal null directions of electro¬ 


day tensor are coincident, just like simple plane waves in 
flat spacetime. These solutions can thus be seen locally 
as generalized plane waves (see the discussion following 
Eq. (22) of Ref. [T^] for more explicit local similarities 
with plane waves), whose propagation directions follow 
ingoing or outgoing shear-free null congruences |107] . In 
a curved spacetime, this implies that they must evade 
being backscattered by the spacetime curvature, or else 
they cannot remain purely ingoing or outgoing. As a con¬ 
crete example, the solution given by Eqs. (22|-(28) has 


only the ingoing wave component cj)Q, while the outgo¬ 
ing wave component (j )2 and the Goulomb background (pi 
vanish identically throughout space and time. To under¬ 
stand how backscattering is avoided, we recall that the 
ingoing solution as specified by (pQ follows the Kinnersley 
tetrad n basis direction, which is tangential to a geodetic 


shear-free null congruence |I08j . By the Goldberg-Sachs 
theorem [98] : 

A strictly nonflat vacuum metric has a multiple principal 
null direction 1°“ iff 1°“ is geodetic and shear free, 
this n direction must also be the direction of two or more 
degenerate GPNDs. We now show explicitly that follow¬ 
ing degenerate GPNDs is responsible for the simplifica¬ 
tions to scattering by spacetime curvature. 

We begin by recalling how electromagnetic waves (with 
allowance for current and charge, so FFE waves are in¬ 
cluded as a special case) scatter off of spacetime curva¬ 
ture. The wave equation satisfied by the Faraday tensor 
is given by [109] 

+Vf,Ja— VaJ&, (49) 


where V'^Vc is the generalized covariant Laplacian oper¬ 
ator. The scattering by spacetime curvature is described 
by the first three terms on the right-hand side of Eq. (491, 


and is a consequence of the tensorial nature of F that al¬ 
lows it to couple to the spacetime curvature through the 
Ricci identities, which when applied to our case gives 


SVlaVtijEcd = RabceF^^d + RabdeFf 


(50) 


and subsequently yields the aforementioned terms. These 
scattering terms imply that, generically, the electromag¬ 
netic waves can propagate inside as well as on the future 
null cone of a light source, as secondary ingoing waves 
can be created by scattering, and so they do not sat¬ 
isfy Huygens’s principle |1()9H ]T2 ] . The remaining terms 
on the right-hand side of Eq. (491 describe scattering by 


charge and current, and are not the scattering we are in¬ 
terested in here. In other words, our consideration in this 
section concentrates on the scattering shared by vacuum 
(without current), FFE (with current), as well as other 
electromagnetic solutions with more generic currents, so 
that our conclusions are not confined to the FFE case. 


magnetic solutions. 
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Since we do not include the stress-energy tensor of the 
electromagnetic field or the plasma in the gravitational 
sector, as per the simplifying convention in FFE compu¬ 
tations we have Rab = 0, and Rabcd reduces to the 
Weyl tensor Cabcd- Both Cabcd and Fab can be written in 
the spinor formalism as m 

Fab = 4>ABeA'B' + eAB^A'B' , (51) 

Cabcd = '^ABCOeA'B'^C'D' +'^A'B'C'D'eABiCD{52) 


Note that as per convention, we have left out the solder¬ 
ing forms like for brevity, with the understanding 

that spinor index pairs like AA' correspond to tensor in¬ 
dex a (same letter but lowercase). The multi-indexed 
spinors can be further factorized into products of their 
respective principal spinors (relating to the principal null 
directions of the original tensors via Eq. E5): 


where ipo is extracted under any dyad ( corre sponding to 
a Newman-Penrose null tetrad via Eq. Ell) that has 0 
as a member of its basis. We will denote such a dyad as 
(0^,C'^). 

We now consider a null electromagnetic wave travelling 
in a purely radiative (Petrov type N) |103[ 1104) space- 
times® with all four GPNDs coinciding. If furthermore, 
the electromagnetic wave travels in the same direction 
as the gravitational wave, so that the doubly degenerate 
EPNDs coincide with the four-fold degenerate GPNDs, 
then all of the contractions in Eq. (56) will vanish, and 
there will not be any curvature scattering term left. Such 
a complete disappearance of scattering can also happen 
when the spacetime is of type III. 

For the Petrov type D Kerr spacetime, the scattering 
term in Eq. (56) does not vanish completely, so there is 
still some residual scattering, but of a simplified structure 


4'AB = 0{Al'B) , (53) 

"^ABCD = 

where more specialized Petrov classes of spacetimes have 
more of the a^s (corresponding to the GPNDs) being 
coincident. 

Now a simple calculation shows that the scattering 
term translates into the spinor language as 

CacbdF'^'^ = ABCd4'^^) ^A'B’ + c.c., (55) 


where c.c. stands for complex conjugation. Substituting 
in Eqs. (53) and (54), we find that the spinor counterpart 


to the scattering term is 


1 


'^abcdcIF'^ = «b) (56) 


where (j, j) are an unordered and unequal pair of num¬ 
bers from {1, 2, 3,4}, while (fc, 1) are an ordered pair con¬ 
sisting of the remaining two numbers, with k > 1. From 
Eq. (56), and recalling that the contraction of coincident 


spinors vanish, it is clear that the more pairs of EPNDs 
and GPNDs that are coincident, the more terms in the 
sum will vanish, leaving us with a scattering term that 
is simpler in its composition, meaning it has fewer inde¬ 
pendent components. 

Such an effect is particularly strong with null electro¬ 
magnetic waves, defined by the property that the two 
EPNDs are coincident or 0 oc i (Sec. 5.1 of [TT3]). This 
implies that 4’ab<I>^^ = 0 due to Eq. (E2). Through 


Eq. (51), this further implies that the two real invariants 


of the Faraday tensor must vanish [57], i.e. Eqs. (23) and 
(24) must be satisfied. In particular, these conditions are 


consistent with the force-free constraints, so there can be 


FFE null solutions. More specifically, using Eqs. (51), 


(53), (Ell) and (E14), it is easy to verify that for a purely 


ingoing null solution, we can write 


Fab = 4>O0A0B^A'B' + c.c. , 


(57) 


CacbdF'^'^ = (j)o^20A0Be-A'B' + C.C. , (58) 

that does not necessarily lead to 6acfcscattering. In par¬ 
ticular, the scattering term only contains 0a and no C,a 
in its spinor form, which helps prevent a contamination 
of Fab by 01 and 02 , as the spinor counterpart of Fab is 
given by 

4>AB = 0O0A0B — 20i0(a0b) + 02CaCb , (59) 

so that 01 and 02 need to multiply with C^a in order to 
pick up spinor and subsequently tensor indices. A rigor¬ 
ous proof for the existence of backscatterless solutions in 
type-D spacetimes is provided by the the general theorem 
(see Refs. UnZlIIIl] and (7.3.14) of Ref. [57]): 

If 1°' is a geodetic and shear-free null congruence and 
analytic, then there is a nonzero solution of the vacuum 
Maxwell equations which is null along C. 

This result applies in the vacuum case and, to a more 
limited extent (restricted to Kerr spacetimes [12]), in the 
force-free case. The proof of the theorem above and the 
discovery of the FFE solutions are rather technical in na¬ 
ture, but we hope that our discussion regarding the con¬ 
traction/annihilation between EPNDs and GPNDs and 
the subsequent simplification/elimination of the scatter¬ 
ing term will serve to help build intuition for these im¬ 
portant results. 


VI. CONCLUSION 

In this paper, we have introduced a new pseudospec- 
tral fully 3D curved spacetime FFE code, with an initial 


Petrov type I: four different GPNDs; type II: two degenerate 
GPNDs and two nondegenerate GPNDs; type D: two sets of 
doubly degenerate GPNDs; type III: a triply degenerate set of 
GPNDs and a nondegenerate one; type N: fourfold degeneracy 
in GPNDs; type O: Weyl curvature tensor vanishes. 



























23 


data solver and an improved constraint damping mecha¬ 
nism. Using this code, we have shown through numerical 
experiments that the backscatterless analytical FFE so¬ 
lutions found by Ref. [T^] are stable against a variety of 
perturbing scenarios, which we selected to avoid restrict¬ 
ing ourselves to special subspaces of the FFE solution 
space. However, with any simulation, one can only con¬ 
strain the growth rate of unstable modes, as the simu¬ 
lations can not be performed for an infinite amount of 
time. We have carried out our simulations for around 
ten light crossing times, and observed the perturbations 
to exit the computational domain after two or three. If 
there exist unstable modes, and they are excited to an 
appreciable initial amplitude, we would expect them to 
remain in the computational domain and then gradually 
grow. This does not appear to be the case. Nevertheless, 
we should be cautious and only place an upper bound 
on the growth rate of any unstable modes assuming the 
worst scenario. Assuming that the unstable modes ex¬ 
ist but are not excited by our initial perturbation at all, 
and instead started at the floor level of numerical noise 
with and AB^ ~ 10“^° (see Fig.j^, our final values 
for these quantities are around 10“^^(see Figs. 11 16 


and 19) after 2000M of simulation, which gives an up¬ 
per bound on the growth rate for VAE'^ and \/ AB"^ at 
around 0.005. For comparison, the fundamental I = 1 
quasinormal mode decay rate for the magnetosphere in 
our study is 0.09. 


In order to carry out concrete numerical studies, we 
had to select particular solutions from the family of in¬ 
finitely many null’*' solutions that can be obtained. We 
have chosen arbitrary representative solutions that share 
all the important physical features listed in Sec. |II C| with 
the rest of their siblings that are not explicitly simulated, 
and these solutions do not have any special features or 
symmetries that would make them more (or less) sta¬ 
ble than the rest. We have also examined both classes 
- the time-dependent “waves” and the time-independent 
“winds” - of the solution family. Therefore, we expect 
the stability conclusion to be generalizable to the ma¬ 
jority, if not all, of the solutions in the parent family. 
In addition, we have chosen the magnitudes of the un¬ 
perturbed solution in such a manner as to ensure the 
tests are carried out within the nonlinear regime of the 
FFE equations, so our results suggest full nonlinear sta¬ 
bility. Furthermore, even though we cannot make math¬ 
ematically rigorous statements as to the asymptotic sta¬ 
bility of the null"*" solutions due to the presence of nu¬ 
merical noise, we note that with all three different types 
of perturbations that we have examined, the perturba¬ 
tions are seen to exit the computational domain after 
two light-crossing times, and the final AB^ and AE'^ are 
small. From a practical point of view, even if the fi¬ 
nal “steady-state” solutions are not exactly the same as 
the null'*' solutions, they would be well approximated by 
them, so that the null'*" solutions can be considered ef¬ 
fectively asymptotically stable. Therefore, despite their 
physical specialness, these solutions are not fragile, and 


can in fact describe physically realistic scenarios like the 
outer magnetospheres of pulsars. We note that one fea¬ 
ture of the FFE null+ solutions that made our numerical 
study possible is the fact that they are globally regular, 
and are therefore amenable to numerical simulations. In 
contrast, the vacuum null solutions in curved spacetimes 
would contain singularities if extended globally [Il[71]. 

Unlike the jet stability (despite many plasma instabil¬ 
ities that can be present in that case, the observed jets 
from active galactic nuclei can extend to several hundred 
kpc), the stability of the null+ solutions is perhaps less 
surprising. The null"'" solutions are more or less isotropic, 
so there is no need to maintain a narrow collimated en¬ 
ergy flux by e.g. ambient gas pressure, and no sharp 
boundaries in the form of a jet surface exist. Therefore, 
we do not have problems like dangerous surface Kelvin- 
Helmholtz modes associated with the vortex sheet on the 
jet surface [sniisi]. Nevertheless, there are many interest¬ 
ing physical properties of the null’*' solutions, such as hav¬ 
ing a null current and being scatter avoiding, which we 
did not know would be stablizing or destablizing. Given 
the stability result from our numerical experiment, one 
can now ask interesting questions such as whether the 
current being null will actually help prevent the onset of 
current-driven instabilities? Would an electromagnetic 
wave propagating close to a degenerate GPND direction 
in fact tend to end up moving exactly along it during 
its journey (in which case the null"'' solutions will really 
have a preferential status)? These will become interest¬ 
ing topics for future studies. 

We caution however, that the stability of the null’*' 
solutions does not necessarily translate into that of the 
null solutions. This is because first of all, the null ex¬ 
act solutions will only satisfy the magnetic dominance 
condition marginally, so any perturbation would likely 
introduce non-FFE influences such as current sheets. In 
the study of jet stability, the presence of current sheets 
at the jet surfaces is long known to encourage instabil¬ 
ities EH, and in a more recent study, it has been seen 
that dynamical current sheets have the tendency to dras¬ 
tically rearrange the underlying force-free solutions even 
when the initial data provided are already a solution to 
the FFE equations Such current-sheet-induced in¬ 

stabilities may well also beset the null solution stability. 
Secondly, taking the limit where the null waves become 
Alfven waves, the background monopole B field in a null"'" 
solution provides a preferred direction for the waves to 
travel in. It is not clear whether the null solutions would 
be more susceptible to changing propagation directions 
when this guidance is taken away. Lastly, it has been ob¬ 
served in the study of jet stabilities that increased mag¬ 
netization has a stabling effect [57!, which can also be 
gone in the null case. As in astrophysical situations, we 
expect a split monopole/dipole-like background field to 
be present, the stability of the null+ solutions observed 
in this work would hopefully mean that they can serve 
as efficient channels for carrying energy across magneto¬ 
spheres unhindered. However, because the background 
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field drops off faster than the null wave component in a 
null'’' solution, eventually, this stability may be lost (if 
the edge of the magnetosphere is not reached first), and 
the ensuing instability may play a role in how the trans¬ 
ported energy finally turns into observational signals. 

Finally, on the analytical front, we have carried out an 
explicit analysis of the scattering of an electromagnetic 
wave by the spacetime curvature, with emphasis on the 
role played by the GPNDs of the Weyl curvature tensor. 
We showed that waves propagating along the degenerate 
GPNDs experience simpler forms of curvature scattering, 
thereby providing some intuition into the perplexing ex¬ 
istence of backscatterless null solutions (although general 
theorems on this subject already exist, their proofs are 
highly technical). One interesting new conclusion is that 
the scattering can vanish completely in a Petrov type 
III or type N spacetime. This would remove a potential 
complication (though may or may not be significant in 
the first place) associated with analysing the coinciding 
electromagnetic counterpart to a gravitational wave, as 
the gravitational wave would not in fact try to scatter its 
electromagnetic companion. The understanding of curva¬ 
ture scattering gained here should also prove useful when 
constructing new analytical solutions to the FFE equa¬ 
tions and other field theories where a lack of scattering is 
desirable, perhaps for the sake of reducing computational 
complexity. 
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(characteristic speeds) in curved spacetime are simply 
given by 


We thank Zachariah Etienne and Matthew Duez for 
useful discussions, and the anonymous referee for many 
helpful comments. The numerical simulations presented 
in this work were performed on the WVU cluster Spruce 
Knob and the Galtech clusters SHG and Zwicky. Spruce 
Knob is supported in part by National Science Eounda- 
tion EPSCoR Research Infrastructure Improvement Go- 
operative Agreement No. 1003907, the state of West 
Virginia (WVEPSGoR via the Higher Education Policy 
Gommission) and WVU. Zwicky is funded by NSF MRI- 
R2 Award No. PHY-PHY-0960291. F. Z. is supported in 
part by National Natural Science Foundation of Ghina 
Grant No. 11443008, Fundamental Research Funds for 
the Gentral Universities Grant No. 2015KJJGB06, and a 
Returned Overseas Chinese Scholars Foundation grant. 


Appendix A: Hyperbolicity of the evolution system 

The study of the characteristic structure of a set of 
evolution equations arises from the need for the initial 
value problem to be well posed. For our pseudospec- 
tral implementation, there is an added urgency because 
the boundary conditions for the overall computational 
domain | 116 H 118 j as well as between the adjacent subdo¬ 
mains [55H5H] are imposed on the characteristic modes. 
Recall that the evolution system can be written as a col¬ 


i>a = Ni/a — /3 • n, (A3) 


while the eigenvector (characteristic mode) expressions 
are unchanged from their flat spacetime counterparts. 
The analysis of the characteristic structure is then es¬ 
sentially independent of the spacetime curvature. 

For the minimal evolution system, the right eigenvalues 
satisfying the equation 


3 ~ CSL 

36 ^ = 


(A4) 


are [75] 


where 




UJ = 


R2 


-N -13-n, 

(AS) 

N — (3 ■ n, 

(A6) 

N{y — w) — /3 • ft, 

(A7) 

N[v -1- w) — /3 • ft, 

(AS) 

-/3 ft, 

(A9) 

-/3 ft, 

(AlO) 

(E X B) 

B2 

(All) 

V(ft-B)2(R2_^2)_ 

(A12) 
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The right eigenvectors are 


ef = (-PE + nxB, PB + fi x e) , (A13) 

e| = (-PE-nxB, PB - ri x E^ , (A14) 

65,4 = (-PB + 1^3 4 n X E + (1 - t'l 4 )B, 

-PE - 4 n X b) , (A15) 

ef = ( 0 , n) , (A16) 

el = ((n • B)n, Pe) , (A17) 

where ^ = 1 / ± to, and for a vector A, (PA) is defined 
as A — (n • A)n. The e“ and e| are fast modes travelling 
at the speed of light. The terms e? and e? are the Alfven 
modes, while e? and e? are the unphysical modes (there 
can only be four physical modes as constraints V • B = 0 
and E • B = 0 reduce the number of independent degrees 
of freedom to four). 

Note that the characteristic speeds become complex 
when > B^, and the evolution system will not be 
strongly hyperbolic. This is one issue that stems from the 
physical constraint of subluminal motion for the plasma 
particles. There is however another hyperbolicity related 
problem with the FFE evolution equations. Namely, even 
when all the constraints are satisfied, we do not have a 
complete set of eigenvectors when (n • B)^ = (n x E)^ 
m- For example, when n • B = 0 and E = Eh, we can 
choose coordinates such that n = x and B = By-, then 
for the minimal system, we have 


n A = 


/O 0 
0 0 
^ 0 
0 0 
0 0 
\0 1 


0 0 

0 0 

0 0 

0 0 

-1 0 

0 0 


0 0 \ 
0 0 
-1 0 
0 0 
0 0 
0 0 / 


(A18) 


whose characteristic equation is 

det(n • A — Al) = A'^(A^ — 1) = 0, (A-19) 


so that we have four zero eigenvalues. In order to allow 
for four eigenvectors corresponding to these zero eigenval¬ 
ues, matrix n • A must have rank 6 — 4 = 2. However, its 
actual rank is 3, so we don’t have a complete set of eigen¬ 
vectors. As an aside, we mention that in some numerical 
implementations, a divergence cleaning scalar field is in¬ 
troduced into the evolution system [53 112011121 ] , which 
enlarges the space of evolved variables, and changes the 
characteristic structure of the evolution system. How¬ 
ever, doing so does not cure this particular hyperbolicity 
problem (details are provided in Appendix [^. Neverthe¬ 
less, the directions for which we do not have a complete 
set of eigenvectors is a set of measure zero among all 
possible fi directions, and there are enough eigenvectors 


to represent the constraint-satisfying solutions even for 
these directions m- However, for constraint violating 
solutions, the constraints may grow on arbitrarily short 
time scales (beyond the ability of our constraint damp¬ 
ing additions to control) when the evolution system is 
not well posed. Thankfully, for the numerical studies in 
this paper, we do not encounter such a situation. 

We note, nevertheless, that it is possible to obtain 
strictly strongly-hyperbolic evolution equations by aug¬ 
menting them with terms that vanish for constraint- 
satisfying solutions, so that only unphysical modes are 
altered. Ref. m provides one such system, and we fur¬ 
ther improve upon it by bringing in more augmentation 
terms and proposing two systems with additional desir¬ 
able properties. In particular, one system remains sym¬ 
metric hyperbolic even when the constraint E • B = 0 
is violated. The other system, although no longer sym¬ 
metric hyperbolic when the constraints are violated, has 
a strongly hyperbolic set of constraint evolution equa¬ 
tions. The details of these augmented systems are given 
in Appendix [B} 


Appendix B: Well-posed FFE evolution systems 

It is desirable for an evolution system to be strongly 
hyperbolic, as then it will be well posed. It has been 
shown in a recent paper m that it is possible to aug¬ 
ment the FFE evolution equations with terms containing 
the derivatives of the constraints, such that the result¬ 
ing system is symmetric hyperbolic when the constraints 
are satisfied. We show in this appendix section that by 
considering additional augmentation terms, it is possi¬ 
ble to make the evolution system retain its symmetric 
hyperbolicity even when the FFE constaint E • B = 0 
is violated. As constraints are never exactly satisfied 
in numerical simulations, such nice off shell (off of the 
constraint surface) properties are obviously desirable. In 
addition, we also provide an alternative augmented evo¬ 
lution system, whose evolution equations do not remain 
symmetric hyperbolic off shell, but whose associated con¬ 
straint evolution equations are strongly hyperbolic and 
particularly simple, so that the constraints evolve in a 
well-understood and controlled manner. 


1. The main evolution equations 

The unaugmented evolution equations are not strictly 
hyperbolic even when the constraints are satisfied. 
Namely when (fi • B)^ = (n x E)^, the matrix hiA'f does 
not possess a complete set of eigenvectors m- This prob¬ 
lem can be cured by adding constraints to the evolution 
equations. Such additions will not change the physical, 
constraint-satisfying solutions, but will modify the char¬ 
acteristic structure of Eq. if the new terms contain 
derivatives. For our case, we consider seven possible ad¬ 
ditional terms that look similar to already existing ones. 
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With them, the evolution equations become 


{dt - C0)E 


{dt - Cp)B 


NK^ + V X (iVB) - (B-VxB-E-VxE - 2KijE^B^ + 2K'E • B + 5E • B) 

^ 1""^ |h j ^ 

- -^TVV • E - aiN— X V(E • B) - a 2 N x E - a 3 iV(E • B)E xV — 

JD^ 

^ 1“"^ 1^ |h ^ 

NK'B - V X (7VE) - aiN -• B - a^N — x V(E • B) - x (7VB) 

B^ B‘^ 

-a^N^E ■ B)B X 


(Bl) 


(B2) 


where we have included the constraint damping term as 
well (it does nothing to the characteristic structure of the 
equations as it does not contain any derivatives). Note 
that in order to acquire additional desirable properties, 
we have considered a larger collection of possible augmen¬ 
tation terms than Ref. m, which included those terms 
whose coefficients are ai, and as, and with these coef¬ 
ficients fixed to 0, 1 and 1, respectively. In other words, 
the minimally augmented system (abbreviated to AU) 
introduced in m corresponds to 

AU : a = (0,0,0,1,1,0,0). (B3) 


in our notation. 

We will propose two additional augmented systems 
AU2 and AU3 given by 


AU2 : a = (-1, 0, -1/2,1,1,0,1/2) 


(B4) 


and 


AU3 : a = (0,1,0,1,1,-1,1) 


(B5) 


respectively. Both of these systems are symmetric hy¬ 
perbolic (see Sec. B 2 of this appendix) just like the AU 
system, but possess additional desirable properties. The 
AU2 system retains its symmetric hyperbolicity when 
E B ^ 0 (see Sec.|B 3[) and the AU3 system has a partic¬ 


ularly simple constraint evolution system (see Sec. B4). 


2. Hyperbolicity of the main evolution equations 
when constraints are satisfied 

The requirement of hyperbolicity of Eqs. |B1| and |B2| 
will imply restrictions on the coefficients a.. To inves¬ 
tigate these restrictions, we first consider the case of 
n • B = 0 and E = Ufi, which is a special case of the 
(n • B)^ = (n X E)^ configurations. Because the curva¬ 
ture of the spacetime impacts the characteristic structure 
of the FEE equations trivially, we will use only flat space- 
time expressions for the rest of this section, with the un¬ 
derstanding that curved spacetime counterparts can be 


recovered using Eq. (A3). For the augmented systems, 
we have that 


n-A = 


0 

0 

E/B 

0 

0 

0 


0 

0 

0 

0 

0 

1 - as 


0 

0 

0 

0 

-1 

0 


0 0 

0 0 

-1 0 

0 


(a4 — a^)E/B 
whose characteristic equation is 

det(n • A — Al) = A^(A^ — 1) = 0, 


(P6) 


(B7) 


so we have four zero eigenvalues. In order to allow for 
four eigenvectors corresponding to this zero eigenvalue, 
matrix n-A must have rank two, which implies 1 —as = 0 
and a 4 — as = 0, or a 4 = 1 = as. We note that all three 


augmented systems given by Eqs. (B3), (B4| and (B5) 
satisfy this requirement. 

We now turn to prove symmetric hyperbolicity of AU2 
and AU3 for generic cases by explicitly calculating the 
symmetrizer S for them. We do so by writing down the 
most general symmetrizer when the constraint E • B = 0 
is satisfied, with each term multiplied by a yet-to-be- 
determined coefficient. We then solve for these coeffi¬ 
cients by ensuring 5'^a(n • A)'>'^ is symmetric for all n. 
This is a tedious but straightforward process. The condi¬ 
tion of a 4 = 1 = a 5 turns out to be necessary to ensure 
the positive definiteness of S^a- The symmetrizer for 
AU2 is simply 


fB^9ij + (C ^)BiBj EiBj + C,BiEj 

^ -B,E, + QE.Bj B^g,, + (C - l)E,E,) 


where C is a free constant, and we require C > 0 for the 
positive definiteness of S. The symmetrizer for AU3 is 


S = 


_ / Ay,, -f (^A - 2A)^ -A^ + 


-A 


B2 




EiBj 

B2 


B2 
Agij H 




B2 
EiEj 
B2 


(P9) 


where A = 1 — U^/U^, and we should pick a ^ > 1/A to 
ensure positive definiteness of S. 
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3. Hyperbolicity of the main evolution equations 
when E ■ B 7 ^ 0 


The a. vector for AU2 in Eq. ( |B4[ ) is chosen to en¬ 
sure that the symmetrizer remains valid when E • B 7 ^ 0 
(this property is not shared by the AU3 or AU evolution 
systems). Indeed, by picking (^ = 1/2, it is straightfor¬ 
ward to verify explicitly that a (n-A)T'^ is symmetric. 
Namely, if we break the greek indices into a pair of spatial 
indices and write a (n • A)^ ^ in a block form 


f {SAEEVk iSAEByk\ 
\{SABE)^k (SABBykJ 


(BIO) 


we then have 


{SAeeW^^^ = 0 = {SABBUe ^^'^, 

{SAEB)ik — {SABE)ki = 0 , 

regardless of the value of E • B. 

In greater detail, the block form of n • A is 

/ {AEBYk iAEByk\ 
yiABEYk {ABBykj 

where for generic a. choices 


(Bll) 


(B12) 


{Aee)" k 


{AEBy k 


{ABEy k 


1, — 


{ABBy k 


1 , — 




(E X B) 


-hfe -I- ai 


(E X n)^' 


Bk+a2- 


E B 






-e^'kfii + X n)fc Ek - 2a3^^(E x nYBk 

, (B X n)-’’ 

kUi + as-^—7^;:^- B 


B 2 


k 5 


(ExB)J' (B X n)J E-B . E • B 

a 4 -- nk + as - — - Ek -I- ag kUi - 2 aj (B x ny Bk ■ 


^2 


B 2 


B 2 


B^ 


(B13) 

(B14) 

(B15) 

(B16) 


Multiplying with symmetrizer S as given by Eq. (B 8 ) with </ = 1/2 and an extra overall factor of 2 for convenience, 
we obtain the components in Eq. (B10[): 


{SAEE)rk = -Bi{E X n)fc -f 2(E x B)jhfc - 2(E x x^iBk + B • (E x n)BiBk + Bi(E x ri)^ 
-2E,(B X h)k + '^^■ B,Bk , 

{SABB)tk = -Ei{B X n)fc -F 2B^{'E x n)^ -f Ej(B x n)fc - 2^j^Bi{B x n)fc - ^ EiEk 


B 2 


+ ^^B • (E X n)EiBk + 2(E x B)^nk + 2(B x h)^Ek - ^ EiEk 

„E-B_ ,, ^ E-BE-(Bxn)^^ 

—2—^^(B X n)iBfc -I ^2 EiBk , 

E B 

{SAEB)tk = 2B'^eikin^ - 2(E x h),Ek + 2-^(E x n)iBk + 2Bi(B x n)fc , 

E • B 

{SABE)ik = —2B^eikin^ — 2(E x h)kEi + 2 (E x n)kBi -|- 2Bk{B x n)i, 


(B17) 


(B18) 

(B19) 

(B20) 

(B21) 


where we have specialized to the AU2 system of Eq. (B4|. We then have {SAEB)ik — {SABE)ki = 0, and it is 
straightforward to show that the antisymmetric part of the diagonal blocks are 


iSAEEhke^^^ = {SABBhke^^^ = 2(B x (E x fi) + n x (B x E) + E x (fi x B))^' = 0 


(B22) 


For completeness, we explicitly write out the charac¬ 
teristic modes and speeds for the AU2 system. The right 
eigenvalues satisfying the equation 


are 


"1 

vs 


in^A'■)'^ = v^el 


(B23) 


- 1 , 

(B24) 

1 , 

(B25) 

V — UJ , 

(B26) 

V + 0J , 

(B27) 

V, 

(B28) 

2v, 

(B29) 
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where v and w are given in Eqs. (All) and (A12|. The 
right eigenvectors are 


Cn = 


et-i = 


3.4 


(-PEH 

(-PE- 

(-PB 


n X B, 
n X B, 


PB + n X E) 
PB — n X E) 


'^3,4 


n X E + (1 - |)B, 


— PE — X B 


ef = (E(B-h)-B(E-n), - E x B:/) , (B33) 

ef = (B,0). (B34) 

We note that the unphysical modes for the AU2 system as 
given by Eqs. (B33) and (B34) are much less complicated 


than for the AU system as given in Ref. m- This is ben¬ 
eficial for inverting the characteristic modes in order to 
obtain the fundamental variables E and B, which is nec¬ 
essary for some pseudospectral implementations such as 
ours, where both the internal (between the adjacent sub- 
domains) [55H55] and the external boundary conditions 
|116H118) are imposed on the characteristic modes, and 
so need to be translated into the fundamental evolution 
variables before they become useful. 

The left eigenvalues satisfying 


ea(.niA ) ^ = 1/ Bp 


(B35) 


are identical to the right eigenvalues in Eq. (B24|-(B29|, 
while the left eigenvectors are 


= 


e„ = 


= (E — n X B, 

—B — n X E 

(E -I- n X B, 


—B -I- fi X E - 

n B 


(R 2 -E^)n-B 
l + L' 


(R2 - e 2) n • B 


1 - 


n 


n 


R2 


E X B 


n B 




R 2 

n B 


B 


— E(a; — -I- fi X B) 


e„ = - 


n B, 


-E X B 


U} 


e„ = 


' n-B 

—E(a; — — fi X B) 

(0,n) , 




^ n B„ 


= (B.E) . 


(B36) 

(B37) 

(B38) 

)■ 

(B39) 

(B40) 

(B41) 


These eigenvectors are degenerate when E = ±n x B, 
in which case we can find alternative complete sets of 


eigenvectors. When E = fi x B, we have 
and we can pick qTn to construct 

= ^3 = 

ef = (4 -fi X q) , 

(B42) 

ef = (n X q, q) , 

(B43) 

ei = (B X (n X q), B x q) , 

(B44) 

el = (-B X 4 B X (n X q)) , 

(B45) 


while the remaining eigenvectors are still valid. When 
E = —fi X B, we have = ^4 = Ij S'lid can use the new 
vectors 


(B30) 

(B31) 

e“ 

= (q, fi X q) , 

(B46) 

^ 4 “ 

= (fi X q, -q) , 

(B47) 



= (-B X (n X q), B x q) , 

(B48) 

(B32) 

et 

= (B X 4 B X (fi X 4) , 

(B49) 


together with the remaining eigenvectors that are still 
valid. 

Lastly, as an aside, we note that for the AU2 system, 
we can further use the identity 

B E B B , 

-xV- =-^xV(E-B) 

B B R2 V y 

-^(E.B)BxV^ (B50) 
to combine terms in the evolution Eqs. (IBII) and 


4. The constraint evolution equations 

We can derive the evolution equations of the con¬ 
straints V • B and E • B from the main evolution Eqs. 
and IB2I The result is 


dtW-B = -a 4 V(vV-B)-a-V(E-B) 

-h4'E-B, (B51) 

at(E-B) = -agV • (v(E • B))-h Oiv • V(E • B) 

-)-4''(E • B) , (B52) 


where 


E X B 

52 


(B53) 


B 1 

(«5 ~ X ^ + (Q^6 + X B(P54) 


4- 

4-' 


-{ae + ar) ( ) • (V x B) 

(05 - 07 - 1 - a 3 )(E X B) • V ^ 


(B55) 


+ (^5 - 02 ) 
-{ae + as) 


B • V X E 

52 

E • V X B 

B 2 


(B56) 


It is desirable for the constraint evolution equations to 
be strongly hyperbolic, so that the constraints evolve in 
a predictable and controlled manner. Such a property 
is especially useful when the main evolution equations 
are not symmetric hyperbolic off shell, as a well posed 
constraint evolution system would signal that at least 
some good behaviors are retained off shell. The n • A 
matrix for the constraint evolution system is simply 


a 4 V ■ n a ■ n 

0 (a5-ai)v-n: 


(B57) 
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which does not have a complete set of eigenvectors when 
V ■ n = 0 unless a • fi = 0, which can be achieved by 

setting = —ag = as in the AU3 system. Note that 

the AU and AU2 systems do not satisfy this condition, 
so their constraint evolution equations are not strongly 
hyperbolic. 

The 9t(E • B) equation simplifies further when ai = 0 
and dr' = 0. So under the coefficient choice = —ag = 
0 ^ 7 , Oi 2 = Q! 5 , and ai = 0 = ag, the evolution equations 
reduce to a pair of decoupled advection equations 

5t(V-B) = -V(a 4 vV-B), (B58) 

at(E-B) = -V(a 5 vE-B). (B59) 


When combined with the = 1 condition for 

the hyperbolicity of the main evolution equations, we 
obtain the a coefficients for the AU3 system as given 
by Eq. ([Ml. 

We note that Eq. (B51) contains a second derivative of 
B whenever 04 ^ 0. Since hyperbolicity requires a 4 = 1, 
this is always the case (shared by all of the AU, AU2 
and AU3 systems ), with both positive and negative im¬ 
plications: A disadvantage of this term is that the sec¬ 
ond derivative increases the sensitivity to high-frequency 
noise in B. However, the Q! 4 -term in Eq. (B51) will cause 
any constraint violations V • B ^ 0 to propagate along 
V, thus allowing it to propagate away. 


that lead to Eq. (AI 8 I, we have a 7 x 7 matrix whose 
characteristic equation is A^(A^ — 1 )^ = 0 , so there are 
three vanishing characteristic speeds. The matrix needs 
to have a rank of 4 for there to be three independent 
characteristic modes of vanishing speed, but the matrix 
actually has a rank of 5. Therefore, the enlarged evolu¬ 
tion system is still not strictly strongly hyperbolic. 


Appendix D: The null solutions in Kerr-Schild 
coordinates 


In our simulations, we specify the metric using the 
Kerr-Schild slicing and coordinates, and the transforma¬ 
tions between the Kerr-Schild coordinates (t, x, y, z) and 
the ingoing Kerr coordinates (z/, r, 0, ijj) are simply (we 
have kept a in case readers need it elsewhere, for this 
paper, one should set a = 0 in all these expressions) 


t = i> — r, X + iy = {r + ia)e'''^ s\a.9, z = rcos 6 *(Dl) 


Appendix C: Hyperbolicity of an evolution system 
with a divergence cleaning scheme 

In our FEE code implementation, it turns out that 
just as observed in Ref. m, the intrinsic accuracy of the 
pseudospectral code is sufficient to keep the V • B = 0 
constraint under control (provided we remove it from the 
initial data of course), and there is no need for any addi¬ 
tional constraint cleaning procedures. We note, however, 
that Refs. |12011121) used a dynamical divergence clean¬ 
ing scheme. Namely, one adds a —NWcj) term to i9tB, 
where is a scalar field satisfying the evolution equation 

{dt - Cp)(t) = -NV ■ B - Na2(j) • (Cl) 

The (f) field then damps the V • B constraint with a time 
scale of (T 2 1 m- 

The introduction of a new field like (j) enlarges the 
space of evolution variables, and would generally alter 
the characteristic structure of the evolution system. Un¬ 
fortunately this does not remove the hyperbolicity prob¬ 
lem when (fi • B)^ = (n x E)^. The characteristic matrix 
for the now enlarged system is 

/ {AEEy k i^EBYk 0 \ 

h • A = (ABEYk (ABByk W , (C2) 

V 0 fik 0 J 

where the four sub-blocks in the top left are the same 
as the minimal system. Under the same assumptions 


with the inverse transformations being 


\J\ j , 

(D2) 

2 1 2 1 2 2 

X + y + z - a , 

(D3) 

arccos , v = i + r , 

(D4) 

“'=‘“(|)+i(i-o) 

1 — arctan (D5) 


Aside from the null solution, we will also add in a 


monopole contribution according to Eq. (32) to ensure 


magnetic dominance, whereby the Faraday tensor be¬ 


comes 


Fab = {(j)om[anb] + 4>iirn[arnb] + n[alb])) ■ (D 6 ) 


Note that the second term for (j)i in the expression 
above vanishes when is given by Eq. (32), which is 
purely imaginary, and so the above expression reduces to 
Eq. ([33]). 


We also note that the Kinnersley tetrad in the Carte¬ 
sian (vector components are in the Cartesian basis) Kerr- 
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Schild coordinates is given by 


f 

A 

(D7) 

l^ + iP 

= sin0 ^^^(r-1-ia) + 1^ 

, (D8) 

l^ 

= cos 9 , 

(D9) 

n° 

A 

(DIO) 

+ in^ 


(Dll) 


A 

= 

(D12) 

m° 

= —^ia sin 9 . 

(D13) 

m} 

= — (3?A cos 9 — i3A) , 

(D14) 

mP 

= — (3A cos 9 + ORA) , 

(D15) 

mP 

= -^rsin0, 

V2 

(D16) 


where A = {r ia) expliip). 

For use in Sec. |IVB[ we also summarize the relation¬ 
ship between the Boyer-Lindquist (t, r, d, (f>) and Kerr- 
Schild (t,x,y,z) spatial coordinates. The transforma¬ 
tions between them are 


t = t + 2Mln 


r 


X + iy = {r + m)exp(i^) sin 0 
z = rcos9 


(D17) 

(D18) 

(D19) 


where 



(D20) 


which in the a = 0 limit gives the Jacobian 


^xks 

3xbl 


/ ^ 2M Q 0 \ 

0 cos(/)sin0 r cos (j) cos 6 —rsm(j)sm6 
0 sin (j) sin 9 r sin (j) cos 9 r cos (j) sin 9 
\0 cos 9 —rsin6> 0 / 


Appendix E: The spinor formalism 


and primed indices (e.g. G W’), while elements of W 
have unprimed indices and no special overhead symbols 
(e.g. G W). We then map between W {W) and its 
dual space W* {W'*) using an antisymmetric spinor cab 
(and ca'ba where it is customary to leave out the over¬ 
head bar on ca'B')- In other words, we raise and lower 
spinor indices with an antisymmetric second rank spinor 
as 

e^BA=U&W\ (El) 


rather than with a symmetric one. The result is that 
spinor self-contractions vanish, i.e. 

= 0. (E2) 

In fact, we have the stronger result (Proposition (2.5.56) 
in Ref. |113j ) that: 

uaP^ = 0 iff a A and Pa are (complex) scalar multiples 
of each other. 

When this happens, we will call oa and Pa coincident. 

We also have a map between the tangent space of the 
spacetime and the tensor product space W ® W given 
by 




V = V 


AA' 


AA' 


where cr are the soldering forms satisfying 

_AA' b _ c 

^AA' — 


b 

-'a 7 


_a _ S. B’ 

(^AA'^a — —^A ^A' 


(E3) 

(E4) 


where the minus signs are due to our metric signature 

choice of (—-|--|—1-) instead of the customary (H-) for 

spinor calculations. The consequence is that whenever we 
translate a contraction between a pair of spacetime in¬ 
dices into the contraction between the corresponding pair 
of double (one primed and the other unprimed) spinor in¬ 
dices, and vice versa, we should add an extra minus sign. 
This step would not be necessary had we adopted the 

(H-) signature as Refs. [571 111311122] did. We note 

that the spacetime null vectors map to particularly nice 
factorized spinor forms 


ja ^A^A' 

' =^AA'^ ? 


(E5) 


With the relatio nshi p (E5) and after applying Eq. (|E4| , 
we see that Eq. (E2) is equivalent to null vectors having 
zero norms under the Lorentzian metric. 

Using the soldering forms, we can also transfer higher 
rank spacetime tensors into their corresponding multi¬ 
spinors. For example the metric maps to 


^ab _a ^AA' BB' 

9 = ^aa'(^bb'9 


(E6) 


For a comprehensive introduction to spinors, please 
consult, for example, Refs. [97l 111311122j . Here for the 
sake of completeness we give a brief summary. 

When spinor bundles can be defined on a spacetime 
(see Ref. ma), we have a two-dimensional complex vec¬ 
tor space W, as well as its complex conjugation W', over 
each spacetime location. The elements of W are written 
with an overhead bar (signifying complex conjugation) 


where 


gAA'BB' _ _^AB^A'B' 


(E7) 


More importantly for us, the Weyl tensor also has a 
spinor counterpart abcd defined by 

n iT( ^ ^ ^AA'BB'CC'DD' 

Gabcd = ^ABCD^A’B'^C'D'CCa d}, Cr^ ( 7 ^ 

^ ^ „AA'BB'CC'DD'(^o\ 

+ 'i’A'B'C'D'eABiCDCba 
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where '^abcd has a much more straight-forward rela¬ 
tionship with its principal spinors than the original tensor 
version Cabcd had with its principal null vectors. Specif¬ 
ically, the algebraic closedness of the complex numbers 
field underlying spinors ensures that we always have the 
factorization 


'^ABCD = , (E9) 

where the a^-'^s are the principal spinors, and then the 
GPNDs of Cabcd are simply given by 

, («) e {(1). (2), (3), (4)}. (ElO) 

Many other quantities take on more transparent forms 
under the spinor formalism. For example, let a spinor 


dyad (0, t) be defined such that it relates to 
Penrose null tetrad by 

a Newman- 

ja ^A-zA' ,A-A' 

1 — ^ AA'^ 0 5^ — ^AA'^ ^ 5 

^A-rA' -a ,A-zA' 

771 — ^ AA'^ ^ 1 ^ AA'^ ^ 

(Ell) 

Then the definitions for the Newman-Penrose scalars un- 

der that tetrad 


^0 = CabcdCmH^m‘^ , 

(E12a) 

^1 = CabcdCn^Crn‘^ , 

(E12b) 

4'2 = CabcdCml^rfCrT^ , 

(E12c) 

4'3 = CabcdC’n^nfn'^ , 

(E12d) 

4'4 = Cabcdrfm'^rPrh'^ , 

(E12e) 


can be rewritten as 


'ko = ABC00^0^0^0^ , (E13a) 

'ki = abcd0^0^0^, (E13b) 

4^2 = ABCD0^0^1-^ , (E13c) 

'^3 = ^ABC00^1-^1-^ , (E13d) 

'^4 = 'if abcdI'"^!'^ , (E13e) 


whereby the decreasing number of times that 0^ ap¬ 
pears in the definitions establishes the hierarchy of decay 
rates of these scalars under the peeling theorem [EnilST], 
a feature not as visible in the tensorial expressions in 
Eq. (E12). Lastly, to carry out some calculations in the 
main text, we note the fact that under the dyad basis, 
EAB can be written as 


CAB = 0ALB — l'A0B ■ (E14) 
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